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CALCULATIONS  OF  SECOND  MOMENTS  OF  NQR 
SIGNALS  FOR  EXPLOSIVE  MATERIALS 


1.  INTRODUCTION 

Over  the  last  few  decades  numerous  terrorist  acts  involving  explosive  materials  have  been 
executed.  These  present  a  real  threat  for  society.  The  problem  of  security  in  planes,  airports, 
postal  service,  and  other  ground  installations  has  become  such  an  important  issue  to  lead  the 
government  to  sponsor  research  on  detection  of  explosive  materials. 

In  recent  years,  the  problem  of  explosives  detection  has  been  addressed  by  the  research 
group  at  the  Naval  Research  Laboratory.  The  work  is  sponsored  by  the  Federal  Aviation 
Administration  and  the  Department  of  Defense,  Office  of  Special  Technologies.  Techniques 
based  on  Nuclear  Quadrupole  Resonance  (NQR)  spectroscopy  have  been  successfully  developed 
>r  the  detection  of  military  type  explosives,  and  narcotics.  The  research  is  aimed  at  the 
detection  of  small  amounts  of  an  explosive,  possibly  hidden  in  luggage  or  a  postal  package, 
where  the  properties  of  the  explosive  material,  or  its  distribution,  make  X-ray  methods  fail. 

In  this  report  we  present  the  calculations  of  the  second  moments  of  *'‘N  NQR  signals  for 
some  military  explosives.  The  calculations  have  been  performed  for  RDX  and  /3-HMX. 

First  we  describe  the  crystal  structures  of  the  explosive  materials.  These  structures  were 
reported  in  the  literature  and  obtained  from  X-ray  and  neutron-diffiraction  analysis.  In  the 
following  two  sections  we  discuss  the  basic  concepts  of  the  vector  and  tensor  analysis  in 
nonorthonormal  systems,  and  outline  of  the  procedure  for  the  calculations  of  second  moments. 
These  two  sections  contain  the  formulae  necessary  for  the  numerical  computations.  A  short 
introduction  to  the  numerical  calculations  follows.  Next,  the  results  of  the  calculations  of  the 
second  moments  are  presented  and  discussed.  The  last  section  summarizes  the  results  of  the 
work  presented.  Some  mathematical  derivations  are  given  in  the  Appendices  A  and  B.  The 
Appendix  C  provides  listing  of  the  most  important  programs. 

2.  CRYSTAL  STRUCTURE  OF  EXPLOSIVE  MATERIALS 

The  well-known  explosive  compound,  cyclotrimethylene-trinitramine,  (C3H(jN(i06),  is 
commonly  called  RDX.  It  exists  in  two  polymorphic  forms,  however,  the  crystal  structure  of 
only  one  form  is  known.  The  structure  determination  of  the  other  form  is  impossible  because 
of  its  extreme  instability. 

Manuscript  approved  November  20,  1995. 
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Another  explosive  compound,  cyclotetramethylene  tetranitramine,  (C4H8Ng08),  is  known 
as  HMX.  It  is  a  by-product  of  the  manufacture  of  RDX.  Four  polymorphic  forms  of  HMX  exist, 
and  of  these  ^-HMX  is  stable  at  room  temperature.  Both  compounds  form  rings  of  alternate 
CH2  and  N-NO2  groups. 

2.  CRYSTAL  STRUCTURE  OF  RDX 

The  crystal  structure  of  RDX  has  been  determined  by  X-ray  diffraction  methods. 
Subsequently,  it  was  refined  from  single-crystal  neutron-diffraction  data' ^  The  unit  cell  is 
orthorhombic  and  space  group  is  Pbca,  Figure  1.  The  unit  cell  parameters  are:  a  =  13.182(2), 
^  ~  11.574(2),  c  =  10.709(2)A.  There  are  eight  molecules  per  unit  cell.  Carbon  atoms  have 
tetrahedral  symmetry,  the  N-NO2  groups  are  planar.  The  molecule  of  RDX  has  a  plane  of  nearly 
mirror  symmetry.  The  plane  is  defined  by  atoms  C2,  Nl,  and  N4. 


Crystal  Structure  of  /3-HMX 

The  crystal  structure  of  j8-HMX  has  been  determined  by  X-ray  and  neutron 
diffraction^’'‘’^  The  compound  crystallizes  in  the  monoclinic  space  group  P2j/c.  There  are  two 
molecules  per  unit  cell  and  the  molecules  are  centro-symmetric.  Figure  2.  Parameters  of  a  unit 
cell  are:  a  =  6.54,  b  =  11.05,  c  =  8.70A,  b  =  124.3°. 
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Figure  2  The  crystal  structure  of  /S-HMX. 


3.  VECTOR  AND  TENSOR  ANALYSIS  IN  THE  NONORTHONORMAL  SYSTEMS 
The  most  important  concept  of  tensor  analysis  is  that  the  properties  of  any  physical  or 
geometrical  quantity  remain  invariant  under  transformation  of  the  coordinate  system  in  which 
its  properties  have  been  originally  described.  In  particular,  considering  crystal  structures,  the 
distance  between  atoms,  angles  between  bonds,  and  the  volume  of  a  unit  cell  must  be 
unchanged.  When  noncoplanar  vectors  span  a  three-dimensional  space  (i= 1,2,3),  then  a  set 
of  reciprocal  vectors  a,*  exists,  such  that 


There  are  two  different  conventions  for  a  definition  of  crystal  parameters  in  the  crystallographic 
literature*.  In  most  textbooks,  the  direct  unit  cell  parameters  (axes  and  angles  between  them)  are 
denoted  a,  b,  c,  a,  b,  g\  the  reciprocal  unit  cell  parameters  are  denoted  a*,  b*,  c*,  a*,  b*,  g*. 
The  basic  relations  between  the  basis  vectors  are: 
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a 


* 


bxc 
V  ’ 


h 

c 


(2) 


V=a'{bxc),  V*  =  V-^  , 


where  V  is  volume  of  a  unit  cell  in  the  direct  lattice,  and  V*  is  volume  of  a  unit  cell  in  the 
reciprocal  lattice.  A  metric  tensor  for  the  direct  lattice  is 


i  = 


u'a  a‘b  a'C 
b-a  b'b  b'C  , 
C'a  C'b  C'C 


(3) 


where  each  tensor  element  is  a  scalar  product  of  the  corresponding  basis  vectors.  Similarly,  the 
metric  tensor  for  reciprocal  lattice  is 


g 


* 


a*-a*  a*-b* 

b*-a*  b*'b*  b*‘C* 

c*-a*  c*-b*  c*-c* 


(4) 


In  our  calculations  we  use  another  convention’.  First  we  define  a  vector  in  the  direct 
basis  and  in  the  reciprocal  basis  of  the  coordinate  system: 


X  =  , 

X  =  x,a< , 


where  a,  are  base  vectors  and  o'  are  reciprocal  base  vectors,  and  i=l,2,3  for  three-dimensional 
lattices.  Subscripts  denote  covariant  quantities,  superscripts  denote  contravariant  quantities,  and, 
according  to  the  Einstein  notation,  summation  is  implied  over  the  repeated  indices.  The  covariant 
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quantities  transform  as  the  base  vectors,  while  contravariant  quantities  transform  as  the 
reciprocal  base  vectors.  Transformation  between  base  and  reciprocal  base  vectors  is  given  by 

a*=g^aj , 

where  g*,  and  g'*  are  two-covariant  and  two-contravariant  metric  tensors,  respectively. 
Transformation  of  vector  components  is  performed  as  follows 


(7) 


In  the  following  formulae  we  list  the  most  important  relations  between  the  metric  tensors: 


Sit  ~ 


g 


y 

y  = 


a 


g'j  = «/  =  "  ■V  > 

«ji  •  =  <5,*  • 


(8) 


Volume  of  a  unit  cell  in  either  basis  can  be  expressed  in  terms  of  determinant  of  the  appropriate 
metric  tensor: 


Scalar  products  include  a  metric  tensor  and  can  be  calculated  as  follows: 

(XY)  =  =  g’^x^^  =  x’^y„  . 

The  following  equations  define  vector  products: 

[XY]^  =  . 


(10) 


(11) 


where  e  is  Levi-Civita  symbol,  and 
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^rmn 

^rmn 

=  0-V2  . 

o  > 

(a) 

^mn 

* 

II 

^rmn 

1 

II 

(b) 

^rmn 

=  0; 

^rmn 

# 

O 

II 

(c) 

(a)  rmn  cyclic  permutation  of  1  2  3,  (b)  rmn  cyclic  permutation  of  3  2  1,  (c)  for  any  two 
indices  equal.  Attention  should  be  brought  to  the  fact  that  when  vectors  X  and  Y  are  defined  in 
the  direct  base,  then  the  resultant  of  their  vector  product  is  in  the  reciprocal  base,  and  vice 
versa.  If  needed,  the  indices  may  be  raised  or  lowered  by  an  appropriate  metric  tensor  according 
to  equations  (7). 

The  next  problem  is  the  calculation  of  bond  lengths.  Let  AT  be  a  vector  connecting  atoms 
A  and  B: 

(13) 

then  the  coordinates  of  the  vector  X  are  given  as  the  difference  between  fractional  coordinates 
of  atoms  A  and  B: 


.  (14) 

The  square  of  the  length  of  the  vector  X  is: 

=  (15) 

which  is  equivalent  to: 


^11 

S12  Si3 

X 

>-* 

821 

822  823 

831 

832  833 

(16) 


An  angle  between  vectors  X  and  Y,  connecting  atoms  B  and  A  and  atoms  B  and  C,  is 
calculated  according  to  the  following  formulae: 
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(17) 


X  =  * 

Y  =  y^af,  =cy^~By' * 


cos  6 


gjjxy 


Many  physical  properties  are  described  by  second  order  tensors,  in  particular  by 
ellipsoids.  Thus  the  last  topic  we  address  here  is  the  computation  of  eigen-values  and  eigen¬ 
vectors  of  a  quadratic  form.  Let  us  consider  an  eigen-problem  for  a  two-covariant  tensor  /3.  In 
the  oblique  coordinate  systems,  the  Lagrange  multiplier  X  is  multiplied  by  the  corresponding 
metric  tensor,  two-covariant  in  this  case: 


where  ^  are  the  coordinates  of  the  eigen-vectors,  which  will  be  computed  later.  The  above 
equation  may  be  simplified  by  lowering  indices  of  jc: 


X^  =  g”^x„  . 

After  substitution  we  obtain  the  following  equation: 


(19) 


(20) 


For  a  three-dimensional  space  the  equation  (20)  r^resents  a  set  of  three  linear  homogenous 
equations  for  x„.  This  set  has  nontrivial  solutions  only  if  the  determinant  of  the  coefficients 
vanishes: 
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(21) 


=  0 , 

Further,  this  leads  to  a  cubic  equation  in  X.  In  practice,  we  deal  with  physical  quantities  that  are 
expressed  by  positive  deflnite  quadratic  forms.  Thus  the  solutions  of  equation  (21)  are  real  and 
positive.  Substitution  of  X^*^*  into  equation  (20)  gives  a  set  of  x„  that  determines  an  eigen-vector. 
For  an  ellipsoid,  eigen-vectors  point  along  principal  axes,  and  eigen-values  determine  their 
lengths. 

4.  CALCULATIONS  OF  SECOND  MOMENTS  IN  NQR  SPECTRA 

In  this  section  we  outline  the  procedure  for  calculating  second  moments  of  NQR  lines  for 
nuclei  with  spin  7=1,  that  was  published  by  Shimon  Vega*.  The  equations  presented  here  can 
be  found  in  that  reference  along  with  detailed  derivations.  Some  equations  are  corrected  for 
errors.  We  limit  our  discussion  to  a  case  where  nuclei  possess  integer  spins,  asymmetric  electric 
field  gradient,  and  no  external  magnetic  field  is  present.  Thus,  the  energy  of  the  system  is  given 
by  a  total  Hamiltonian  ^—9^  +9^  (i.e.  sum  of  quadrupolar  and  dipolar  Hamiltonians). 

Several  contributions  to  the  second  moment  are  expected.  We  consider  interactions 
between  completely  equivalent  nuclei,  interactions  between  nuclei  with  the  electric  field  gradient 
(EFG)  tensors  of  the  same  magnitude  but  different  orientation,  and  also  we  calculate  second 
order  contribution  from  nuclei  with  spin  7  =  1/2  (*H). 

NQR  spectral  lines  may  be  either  inhomogeneously  or  homogeneously  broadened.  The 
inhomogeneous  broadening  arises  if  there  is  a  distribution  of  quadrupolar  coupling  constants. 
Such  a  distribution  may  occur  when  impurities  are  present  in  a  crystal,  and  also  when  defects 
in  the  crystal  lattice  are  abundant.  Both  factors  will  affect  the  EFG  experienced  by  a  nucleus 
contributing  to  the  NQR  signal;  the  effect  will  be  most  prominent  in  the  ionic  crystals  where 
impurities  and  crystal  imperfections  alter  significantly  the  electronic  structure  of  the  adjacent 
ions. 

In  the  covalent  crystals  these  effects  are  less  significant,  and  the  broadening  origins 
mainly  in  dipolar  interactions  with  neighboring  nuclei  in  a  molecule,  which  is  homogeneous  in 
nature.  Thus  it  is  convenient  to  characterize  the  line-shape  by  second  moments. 
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The  n*  central  moment  of  a  distribution  /((a)  is 


{<0")  =  f<0"/(0))<i(u)  /  //(«)<i(u)  , 

0  a 

and  the  n*  moment  related  to  the  center  of  a  distribution  is 


oo 

/( 


(22) 


(23) 


In  our  case  /(a))  is  the  line-shape  function  of  a  resonance  signal: 

/(«)  =  .  (24) 

a,b 

where  w  is  frequency  in  radians  per  seconds,  7  is  the  gyromagnetic  ratio,  fi  is  Plank‘s  constant, 
Ir  is  in  the  direction  of  rf  radiation,  and  represent  energy  of  levels  a  and  b,  respectively. 
Summation  is  carried  over  all  states  present  in  the  system.  Thus,  in  general,  several 
spectroscopic  signals  are  included,  while  the  shape  function  of  one  line  at  a  time  is  desired.  To 
make  such  a  restriction,  the  Hamiltonian  of  the  system  is  truncated,  which  means  that  the  off- 
diagonal  elements  of  ^  except  those  connecting  the  degenerate  states  of  ^  are  excluded. 
The  quadrupolar  Hamiltonian  9{q  is: 

=  E ^<31  '  E  A [34  -It  *n [lit  - ly] )] 

i  i 

(25) 

A,  =e\Q,/ 41,(21, -l), 

where  rj  is  the  asymmetry  parameter,  A  is  the  quadrupole  coupling  constant,  Q  is  the  electric 
quadrupole  moment,  and  K.  are  the  principal  values  of  the  EFG  tensor.  The  components  of 
angular  momentum  are  chosen  along  the  eigen- vectors  of  the  EFG  tensor. 
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The  dipolar  Hamiltonian  ^  is 


=  E^if  =  »>'E  (YiT/'-s )[/,•(,- (3/ rj)(/,T(,)(//rj,)]  ,  (2® 

i  i>j 


^>j- 
>'*/  = 


connecting  nuclei  i  and  j.  Further 

=  r,Yi/rl 

){XOS^I^I^*YOS^ 

^yi^yj 

*W-3WV» 

"  Wl » 

sme^jCos<t>^  , 

xos^  = 

1-3X|, 

YOS^  = 

1-31^^ 

cosd^j  , 

Z05,= 

1-3Z|. 

(27) 


where  6  denotes  angle  between  vector  and  z  direction,  and  <{>  denotes  angle  between  projection 
of  on  the  xy  plane  and  x  direction. 

For  completely  equivalent  nuclei  the  second  moment  is 


(A6)^)  =  -. 


2(2/+l))>^tr[/J^.„  j-j.t 


E«[ 


2 

(2/7-1)2 


(28) 


where  ’  is  a  truncated  Hamiltonian. 
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First  the  energies  and  eigen-functions  of  ^  single  nucleus  are  calculated: 


E,  =  , 

|;c>  =  (l/,S)(hl)- 

1  -i» . 

1,)  =  (l/^)(hl>  + 

i-i».  <“' 

=  -2A, 

|z>  =  |0) , 

Then  this  basis  set  is  used  for  computation  of  energies  and  eigen-functions  of  and  (for 

two  nuclei  i  and  j) : 

+ 

II 

II 

II 

\yx)=\y\  \x)j , 

\xy)^\x\\y)j , 

II 

1 

\xx)=  |j;>i  \x)j  , 

(30) 

II 

II 

1 

1 

\yz)=\y)i  \z)j , 

\zy)=\z)i\y)j » 

tl 

II 

1 

+ 

\xz)=  \z)j , 

\zx)=\z}i  \x)j , 

II 

|Z2)=  |z),.  \z)j , 

Elements  of  the  truncated  Hamiltonian  are  calculated  in  the  new  basis: 

and  for  the  transition  |y)  ^|z)  there  are  twelve  elements  of  -I-  ).  All  of  them  are  equal 

to  one. 

Finally,  after  evaluation  of  the  commutator  and  traces  in  Equation  (28),  the  second 
moments  for  completely  equivalent  nuclei  are: 
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(A<i)^)  (I)')**  |z))  = 


(32) 


where 


(1/%*)  (SXOSl*4YOSl*4ZOSl-SYOS,jZOS^)  , 

J*i 


E 


N  n 


E  E 


(33) 


Here,  the  summation  is  carried  over  all  equivalent  nuclei  in  a  unit  cell  (n),  and  a  total  of  N 
nuclei  are  included  in  the  calculations.  The  meaning  of  this  sum  is  that  there  is  N-1  interactions 
between  one  nucleus  and  all  others. 

Next  we  include  equivalent  nuclei  with  different  directions  of  the  EFG  tensors  residing 
in  two  sites  I  and  H.  The  sum  in  Equation  (28)  must  be  replaced  by  »  where 

J>i  ijc  k>l 


i  and  j  represent  nuclei  in  site  I,  while  k  and  /  represent  those  in  site  H.  Now  there  are  following 
contributions  to  the  second  moment:  (A(i)^)(I,I) ,  (Ao)^)(n,n),  (A(«)^)(I,n),  where 


(A<»)^)(I,I)  =  (Aa)^)cos^nrV(cos^qr*  + cos^qr“) , 

and  a*  and  a“  are  the  angles  between  direction  of  radiation  and  direction  of  7,  for  nuclei  in 
sites  I  and  II. 

The  calculation  of  the  second  moments  for  nuclei  in  two  different  sites  ( I,  II)  requires 
a  revision  of  the  coordinate  systems.  So  far  all  the  angular  momentum  operators  have  been 
written  in  the  coordinate  system  of  one  site,  which  was  common  for  all  nuclei  under 
consideration  and  was  determined  by  a  set  of  eigen-vectors  of  the  EFG  tensor.  In  order  to  apply 
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the  above  formulae,  the  vector  quantities  for  the  nuclei  in  both  sites  must  be  referred  to  the  same 
coordinate  system.  Thus  the  angular  momentum  operators  in  the  coordinate  system  of  site  I  are 
transformed  to  the  coordinate  system  of  site  H: 


=  Sr'  * 


(37) 


where  R„.  is  the  transformation  matrix  between  the  coordinate  systems  defined  by  the  principal 
axes  of  EFG  tensors  of  the  nuclei  in  site  I  and  site  H.  To  further  facilitate  calculations,  we 
transform  coordinates  fix)m  the  crystallographic  axes  system  to  the  orthonormal  system  of  the 
same  system  as  defined  for  the  angular  momentum  operators.  The  derivation  of  both 
transformation  matrices  are  given  in  Appendices  A  and  B. 

The  nonvanishing  elements  of  the  truncated  Hamiltonian  and  (/„•  +  /,*)’  are 


^23  “  ^32  “  ^^xs^ik^ik* 

hse  =  ^65  =  -  ^^xVik  -  '^KVik » 

*15  “  *51  “  *37  “  *73  “  *69  ~  *96  “  COSflT  , 

*16  ”  *61  ~  *28  ”  *82  ”  *59  ”  *95  ”  COSflf  . 


(38) 


The  contribution  to  the  second  moment  from  nuclei  in  sites  I  and  n  is 


(Aa)2)(I,II)  (|y)-|z))  = 

(39) 

-8  COS  hr*  COS  (cos^hr*+cos^O'^  ^23^8  ^  ♦ 
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Dipole-dipole  interactions  between  nitrogen  nuclei  and  protons  quite  often  contribute  to  the 
second  moment  of  an  NQR  signal.  However,  in  the  absence  of  external  magnetic  field  and  for 
nonzero  asymmetry  parameter,  this  contribution  vanishes  in  the  first  order.  Thus,  only  the 
second  order  contribution  is  considered: 


(b)**  U))  = 


{ izoslzosl 


*  9[zosl{xl*Yl)  *  ZOS^(X^y^)] 

*  i\{xlxl*YlYl*XlYl*YlX^  ). 


(40) 


The  above  equation  ends  the  list  of  the  contributions  that  may  be  important  for  the  case  of  nuclei 
with  spin  I  =  1/2,  nonzero  asymmetry  parameter,  and  when  NQR  signals  are  measured  in  tiie 
absence  of  external  magnetic  field. 


5.  NUMERICAL  CALCULATIONS 

Programs  for  the  calculations  of  second  moments  of  NQR  signals  are  written  in 
MATLAB,  Version  4.0,  released  by  The  Math  Work,  Inc.  in  1993.  The  most  important 
programs  are  listed  in  the  Appendix  C.  The  computations  have  been  performed  on  a  personal 
computer.  Gateway  2000. 

Two  files,  rdx.m  and  hmx_b.m,  include  crystallographic  data  from  literature  and  provide 
input  to  all  other  programs.  The  calculation  of  second  moments  of  NQR  signals  also  requires 
knowledge  of  the  directions  of  the  principal  axes  of  electric  field  gradient  tensors.  For  RDX  and 
/3-HMX  the  EFG  tensors  have  not  been  measured.  However,  the  orientation  of  the  tensors  can 
be  determined  from  the  electronic  structure  of  the  molecules  of  RDX  and  /S-HMX’.  Thus,  Y 
axis  coincides  with  N-N  bond  in  N-NO2  group,  X  axis  is  perpendicular  to  Y,  and  it  is  in  the 
plane  of  N-NO2  group.  Axis  Z  is  perpendicular  to  XY  plane. 
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6.  RESULTS  AND  DISCUSSION 

This  section  presents  numerical  results  of  the  second  moment  calculations  for  ^-HMX 
and  RDX.  In  most  computations  the  nuclei  in  the  translation  limits  t=3  were  included. 
Translation  limits  refer  to  the  number  of  unit  cells  in  each  crystallographic  direction  (positive 
and  negative)  taken  from  an  arbitrary  unit  cell  chosen  as  an  origin.  Thus  (2t  +  1)’  unit  cells 
were  considered.  Table  1  shows  an  example  of  how  fast  convergence  is  achieved.  The  direction 
of  jrf  radiation  is  fixed.  The  value  of  stops  changing  for  t>3.  Thus,  for  calculations 

of  second  moments  in  the  first  order  this  limit  was  used. 


Table  1.  Data  for  jS-HMX.  The  convergence  of  the  second  moments  calculations 
versus  the  translation  limits,  t.  The  angles  between  the  direction  of  ff  radiation 
and  the  X  axis  of  the  EFG  tensors  for  sites  I  and  n  are  given  by  cosines  cosi  and 
cosn,  respectively 


t 

cosI 

cosn 

<a^>W 

<w*> 

0 

0.800 

-0.875 

630 

9934 

3606 

119 

1 

0.800 

-0.875 

3613 

10191 

3858 

133 

2 

0.800 

-0.875 

3708 

10219 

3886 

133 

3 

0.800 

-0.875 

3729 

10224 

3891 

134 

4 

0.800 

-0.875 

3736 

10226 

3893 

134 

5 

0.800 

-0.875 

3738 

10226 

3894 

134 

Data  in  the  two  following  tables  are  sorted  according  to  the  largest  contributions  to  the 
second  moment  (SM)  from  the  dipolar  interactions  between  nuclei  i  and  j.  Only  the  largest 
contributions  are  listed.  Note  the  difference  between  SM  and  <  .  SM  is  the  contribution  to 

the  second  moment  from  one  pair  of  nuclei  only,  while  <w^>  is  obtained  by  averaging  over 
all  atoms.  The  nuclei  i  are  in  an  arbitrarily  chosen  unit  cell  taken  as  the  origin.  The  nuclei)  are 
either  in  the  original  unit  cell  or  in  the  neighboring  unit  cell  with  the  translation  parameters  m, 
n,  and  p.  The  distance  between  nuclei  i  and  j  is  fjj.  The  direction  of  tf  radiation  in  the 
orthonormal  system  is  [1  -2  6]. 
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Table  2.  Data  for  /3-HMX.  Contributions  to  the  second  moments  from 
completely  equivalent  nuclei  N(2). 


• 

1 

- 

• 

J 

m 

n 

p 

ru  [A] 

SM  [radVsec^] 

2 

- 

4 

0 

0 

0 

3.935 

4868.26 

1 

- 

3 

0 

0 

0 

3.935 

3242.58 

2 

- 

4 

1 

0 

1 

4.358 

1025.50 

1 

- 

3 

-1 

0 

-1 

4.358 

604.78 

2 

- 

4 

1 

0 

0 

5.008 

235.80 

1 

- 

3 

-2 

0 

-1 

7.071 

202.15 

2 

- 

4 

2 

0 

1 

7.071 

165.75 

2 

- 

2 

1 

0 

0 

6.540 

163.05 

2 

- 

2 

-1 

0 

0 

6.540 

163.05 

1 

- 

1 

1 

0 

0 

6.540 

163.05 

1 

- 

1 

-1 

0 

0 

6.540 

163.05 

1 

- 

3 

-1 

0 

0 

5.008 

73.79 

2 

- 

4 

-1 

0 

0 

9.562 

32.96 

1 

- 

3 

1 

0 

0 

9.562 

29.16 

2 

- 

2 

0 

0 

1 

8.700 

14.77 

2 

- 

2 

0 

0 

-1 

8.700 

14.77 

1 

- 

1 

0 

0 

1 

8.700 

14.77 

1 

- 

1 

0 

0 

-1 

8.700 

14.77 

2 

- 

2 

-2 

0 

-1 

10.887 

12.65 

1 

- 

1 

2 

0 

1 

10.887 

12.65 

2 

- 

2 

2 

0 

1 

10.887 

12.65 

1 

- 

1 

-2 

0 

-1 

10.887 

12.65 

1 

- 

3 

-1 

-1 

-1 

11.355 

8.85 

1 

- 

3 

-3 

0 

-1 

12.906 

5.48 

2 

- 

4 

1 

-1 

0 

11.620 

5.28 

2 

- 

4 

3 

0 

1 

12.906 

5.00 

1 

- 

3 

-1 

-1 

0 

11.620 

4.97 

2 

- 

4 

1 

1 

1 

12.380 

4.88 

2 

- 

4 

0 

-1 

0 

11.199 

4.65 
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Table  3.  Data  for  j8-HMX.  Contributions  to  the  second  moments  from 
equivalent  nuclei  N(2)  with  the  different  directions  of  the  EFG  tensors. 


• 

1 

- 

j 

m 

n 

p 

ru[A] 

SM  [radVsec^] 

3 

- 

4 

1 

-1 

0 

5.774 

238.78 

1 

- 

2 

-1 

-1 

-1 

5.774 

238.78 

3 

- 

4 

1 

0 

0 

5.774 

238.78 

1 

- 

2 

-1 

0 

-1 

5.774 

238.78 

4 

- 

1 

0 

1 

0 

6.608 

137.91 

2 

- 

3 

0 

0 

1 

6.608 

137.91 

4 

- 

1 

0 

1 

1 

6.608 

114.58 

2 

- 

3 

0 

0 

0 

6.608 

114.58 

2 

- 

3 

0 

1 

0 

7.472 

59.85 

4 

- 

1 

0 

0 

1 

7.472 

59.85 

2 

- 

3 

0 

1 

1 

7.472 

52.29 

4 

- 

1 

0 

0 

0 

7.472 

52.29 

2 

- 

1 

1 

1 

1 

5.774 

44.44 

4 

- 

3 

-1 

0 

0 

5.774 

44.44 

4 

- 

3 

0 

1 

0 

7.470 

26.89 

2 

- 

1 

0 

0 

1 

7.470 

26.89 

4 

- 

1 

-1 

1 

0 

7.374 

20.35 

2 

- 

3 

1 

0 

1 

7.374 

20.35 

2 

- 

1 

0 

1 

1 

7.470 

17.65 

4 

- 

3 

0 

0 

0 

7.470 

17.65 

2 

- 

1 

0 

0 

0 

8.571 

16.16 

4 

- 

3 

0 

1 

1 

8.571 

16.16 

4 

- 

3 

0 

0 

1 

8.571 

15.46 

2 

- 

1 

0 

1 

0 

8.571 

15.46 

3 

- 

2 

-1 

0 

-1 

7.374 

14.15 

1 

- 

4 

1 

-1 

0 

7.374 

14.15 

3 

- 

2 

1 

0 

0 

7.374 

14.15 

1 

- 

4 

-1 

-1 

-1 

7.374 

14.15 

1 

- 

4 

0 

-1 

0 

6.608 

13.65 

3 

- 

2 

0 

0 

-1 

6.608 

13.65 
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Table  4.  Data  for  ^-HMX.  Contribution  to  the  second  moment  in  the  first 
order  for  N(2)  for  different  directions  of  if  radiation;  translation  limits  t  =  3; 
angles  between  direction  of  jjf  radiation  and  X  axis  of  the  EFG  tensors  of  nuclei 
in  sites  I  and  11  are  given  by  cosines  cosi  and  cosll,  respectively;  <  >  is  in 

units  of  [radVsec^]. 


cosI 

cosn 

-0.910 

0.413 

2583 

3874 

5675 

no 

-0.940 

0.179 

2482 

4507 

6602 

117 

-0.864 

-0.213 

2264 

4402 

6449 

115 

-0.961 

0.370 

2562 

4067 

5958 

112 

-0.996 

0.122 

2449 

4602 

6741 

117 

-0.907 

-0.286 

2236 

4249 

6225 

113 

-0.933 

0.275 

2528 

4298 

6296 

115 

-0.965 

0.031 

2403 

4666 

6836 

118 

-0.884 

-0.359 

2203 

4010 

5875 

no 

-0.371 

0.759 

2593 

902 

1321 

69 

-0.447 

0.268 

2617 

3436 

5034 

105 

-0.298 

-0.539 

2164 

1094 

1602 

70 

-0.535 

0.892 

2617 

1233 

1807 

75 

-0.894 

0.201 

2498 

4446 

6513 

116 

-0.365 

-0.797 

2188 

810 

1187 

65 

-0.426 

0.499 

2645 

1971 

2888 

87 

-0.555 

-0.152 

2252 

4344 

6363 

114 

-0.324 

-0.870 

2215 

570 

835 

60 

0.625 

0.545 

2127 

2653 

3887 

93 

0.707 

0.106 

2309 

4569 

6693 

116 

0.530 

-0.559 

2647 

2212 

3241 

90 

0.784 

0.537 

2142 

3180 

4658 

100 

0.970 

-0.036 

2406 

4664 

6833 

118 

0.617 

-0.766 

2642 

1838 

2692 

85 

0.686 

0.298 

2194 

3927 

5753 

109 

0.800 

-0.230 

2525 

4315 

6321 

115 

0.566 

-0.844 

2628 

1449 

2123 

79 
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Table  5.  Data  for  j8-HMX.  Contribution  to  the  second  moment  in  the  first  order 
for  N(3)  for  different  directions  of  rf  radiation;  translation  limits  t  =  3;  angle 
between  direction  of  if  radiation  and  X  axis  of  the  EFG  tensors  of  nuclei  in  sites 
I  and  n  are  given  by  cosines  cosi  and  cosll,  respectively;  <w^>  is  in  units  of 
[radVsec^. 


cosI 

cosn 

<w*> 

-0.910 

0.882 

3731 

-0.940 

0.964 

3731 

-0.864 

0.969 

3728 

-0.961 

0.885 

3730 

-0.996 

0.974 

3731 

-0.907 

0.973 

3730 

-0.933 

0.799 

3727 

-0.965 

0.881 

3729 

-0.884 

0.891 

3731 

-0.371 

0.372 

3731 

-0.447 

0.588 

3718 

-0.298 

0.548 

3673 

-0.535 

0.394 

3715 

-0.894 

0.939 

3730 

-0.365 

0.574 

3697 

-0.426 

0.164 

3612 

-0.555 

0.387 

3709 

-0.324 

0.396 

3724 

0.625 

-0.588 

3730 

0.707 

-0.555 

3721 

0.530 

-0.278 

3667 

0.784 

-0.843 

3730 

0.970 

-0.890 

3729 

0.617 

-0.405 

3701 

0.686 

-0.858 

3722 

0.800 

-0.875 

3729 

0.566 

-0.471 

3725 

<«*> 

11586 

4410 

140 

10945 

4166 

137 

9956 

3789 

132 

12157 

4627 

143 

11485 

4371 

140 

10445 

3976 

135 

12958 

4932 

147 

12246 

4661 

144 

11138 

4239 

138 

11213 

4268 

139 

8231 

3133 

123 

5135 

1954 

104 

14562 

5542 

154 

10684 

4067 

136 

6471 

2463 

112 

19587 

7455 

175 

15124 

5756 

157 

9038 

3440 

127 

11906 

4532 

142 

13899 

5290 

151 

17627 

6709 

167 

10430 

3970 

135 

12197 

4642 

143 

15700 

5975 

159 

8759 

3334 

126 

10224 

3891 

134 

13261 

5047 

148 
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Table  6.  Data  for  RDX.  Contribution  to  the  N(l)  second  moment,  to  first  order, 
for  two  directions  of  if  radiation;  translation  limits  t  =  3;  angle  between  the 
direction  of  ^radiation  and  X  axis  of  the  EFG  tensors  of  nuclei  in  two  different 
sites  I  and  II  are  given  by  cosines  cosi  and  cosll,  respectively;  is  in  units 

of  [rad^/sec^]. 


cosI 

cosn 

0.732 

0.885 

2856 

210 

55 

0.704 

0.838 

2849 

232 

56 

0.625 

0.728 

2846 

257 

56 

0.501 

0.565 

2849 

287 

56 

0.362 

0.389 

2862 

320 

56 

0.865 

0.969 

2858 

212 

55 

0.861 

0.939 

2854 

234 

56 

0.783 

0.823 

2853 

260 

56 

0.631 

0.630 

2860 

291 

56 

0.458 

0.422 

2878 

326 

57 

0.959 

0.998 

2858 

213 

55 

0.972 

0.977 

2855 

235 

56 

0.898 

0.859 

2857 

262 

56 

0.732 

0.653 

2868 

294 

56 

0.536 

0.431 

2891 

330 

57 

0.982 

0.956 

2854 

214 

55 

0.993 

0.924 

2853 

237 

56 

0.922 

0.807 

2857 

264 

56 

0.764 

0.615 

2870 

297 

56 

0.575 

0.409 

2899 

333 

57 

0.945 

0.861 

2848 

214 

55 

0.939 

0.812 

2847 

237 

56 

0.869 

0.701 

2852 

265 

56 

0.736 

0.539 

2868 

298 

56 

0.575 

0.365 

2900 

335 

57 

-0.271 

-0.347 

2845 

314 

56 

-0.091 

-0.186 

2843 

330 

56 

0.140 

0.024 

2927 

420 

58 
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Table  7.  Data  for  RDX.  Contribution  to  the  N(2)  second  moment,  to  first  order, 
for  different  directions  of  rf  radiation;  translation  limits  t  =  3;  angle  between 
direction  of  //"radiation  and  X  axis  of  the  EFG  tensors  of  nuclei  in  sites  two 
different  sites  I  and  11  are  given  by  cosines  cosi  and  cosll,  respectively;  <(j?> 
is  in  units  of  [rad^/sec^. 


cosI 

cosn 

-0.592 

-0.483 

1808 

554 

49 

-0.579 

-0.674 

1810 

437 

47 

-0.524 

-0.846 

1856 

325 

47 

-0.432 

-0.957 

1932 

241 

47 

-0.326 

-0.997 

2037 

192 

47 

-0.395 

-0.419 

1806 

535 

48 

-0.363 

-0.632 

1882 

393 

48 

-0.294 

-0.831 

1998 

273 

48 

-0.198 

-0.952 

2117 

201 

48 

-0.102 

-0.985 

2273 

175 

49 

-0.143 

-0.317 

1940 

463 

49 

-0.076 

-0.529 

2180 

331 

50 

0.011 

-0.734 

2348 

236 

51 

0.097 

-0.864 

2446 

200 

51 

0.161 

-0.905 

2511 

202 

52 

0.119 

-0.192 

2716 

521 

57 

0.217 

-0.376 

2756 

407 

56 

0.314 

-0.562 

2753 

309 

55 

0.382 

-0.695 

2740 

295 

55 

0.412 

-0.758 

2753 

301 

55 

0.341 

-0.071 

2521 

758 

57 

0.449 

-0.220 

2743 

569 

58 

0.542 

-0.375 

2783 

468 

57 

0.596 

-0.503 

2815 

459 

57 

0.607 

-0.585 

2834 

444 

57 

0.491 

0.967 

1909 

223 

46 

0.561 

0.922 

1872 

226 

46 

0.622 

0.820 

1845 

232 

46 

0.648 

0.648 

1879 

256 

46 

0.627 

0.428 

1879 

368 

47 
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Table  8.  Data  for  RDX.  Contribution  to  the  N(3)  second  moment,  to  first  order, 
for  different  directions  of  if  radiation;  translation  limits  t  =  3;  angle  between 
direction  of  j;/"  radiation  and  X  axis  of  the  EFG  tensors  of  nuclei  in  two  different 
sites  sites  I  and  11  are  given  by  cosines  cosi  and  cosll,  respectively;  <(j?>  is  in 
units  of  [rad^/sec^]. 


cosI 

cosn 

<<J> 

0.633 

-0.608 

2745 

366 

56 

0.751 

-0.429 

2783 

446 

57 

0.834 

-0.195 

2807 

526 

58 

0.856 

0.053 

2828 

420 

57 

0.825 

0.266 

2897 

456 

58 

0.460 

-0.727 

2776 

262 

55 

0.578 

-0.549 

2950 

318 

57 

0.669 

-0.290 

3083 

376 

59 

0.697 

-0.005 

3083 

664 

61 

0.670 

0.236 

3047 

559 

60 

0.230 

-0.814 

2773 

174 

54 

0.324 

-0.643 

3094 

209 

57 

0.408 

-0.378 

3388 

365 

61 

0.449 

-0.073 

3331 

651 

63 

0.448 

0.185 

3100 

604 

61 

-0.017 

-0.842 

2919 

138 

55 

0.040 

-0.678 

3093 

177 

57 

0.105 

-0.426 

3324 

316 

60 

0.159 

-0.135 

3548 

561 

64 

0.194 

0.121 

3023 

594 

60 

-0.232 

-0.816 

2783 

157 

54 

-0.203 

-0.658 

2907 

225 

56 

-0.156 

-0.433 

2987 

354 

58 

-0.097 

-0.177 

2973 

480 

59 

-0.040 

0.058 

3497 

518 

63 

-0.915 

-0.280 

2694 

382 

55 

-0.942 

-0.410 

2543 

336 

54 

-0.929 

-0.552 

2289 

268 

51 

-0.851 

-0.678 

2531 

222 

52 

-0.711 

-0.754 

2896 

211 

56 
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In  Table  9  we  list  the  powder  average  of  second  moments  calculated  for  ^-HMX  and 
RDX.  The  directions  of  if  radiation  were  represented  by  vectors  from  the  center  of  a  cube  to 
the  nodes  on  its  grid.  The  grid  on  each  wall  was  6x6. 

Table  9.  Data  for  /3-HMX  and  RDX.  Powder  average  of  second  moments. 


jg-HMX 


N(2) 

2418 

3910 

5727 

no 

N(3) 

3660 

10843 

4127 

133 

RDX 

N(l) 

2906 

244 

56 

N(2) 

2243 

454 

52 

N(3) 

2857 

367 

57 

1/2 


The  contributions  SM  to  the  second  moment  from  completely  equivalent  nuclei  and  from 
equivalent  nuclei  with  the  different  directions  of  the  EFG  tensors  decrease  rapidly  with 
increasing  distance  (Tables  2  and  3).  The  angular  factor  influences  the  value  of  SM  by 
approximately  one  order  of  magnitude  for  atoms  separated  by  the  same  distance  (Table  3).  For 
j8-HMX  the  contributions  from  completely  equivalent  nuclei  are  larger  than  those  with  the 
different  directions  of  the  EFG  tensors.  The  values  of  <  ’  and  <  >  "  "  are  larger  than 

the  value  of  which  also  is  related  to  the  distance,  since  intramolecular  distances  are 

shorter  than  intermolecular  ones,  and  only  completely  equivalent  nuclei  are  in  the  same 
molecule.  Thus  these  atoms  contribute  the  most  to  the  second  moment. 

For  j8-HMX  the  value  of  <  is  larger  for  N(3)  than  for  N(2),  which  is  explained 
by  the  difference  in  the  intramolecular  distances  between  atoms  related  by  center  of  inversion. 
They  are  2.77  and  3.94A  for  N(2)  and  N(3),  respectively. 

The  dependence  between  values  of  <  >  *’*,  and  <  on 

direction  of  if  radiation  is  given  in  Tables  4  and  5  for  ^-HMX,  and  in  Tables  6,  7,  and  8  for 
RDX.  The  angular  factor  changes  the  value  of  by  less  than  two. 

In  comparison  to  jS-HMX  the  second  moments  <(!)*>*'*  for  RDX  are  significantly 
smaller.  Especially,  the  contribution  from  the  completely  equivalent  nuclei  are  almost  irrelevant. 
Again,  this  is  due  to  the  large  distances,  since  there  is  no  equivalent  nuclei  in  the  same 

molecule. 
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The  linewidths  of  ‘^N  NQR  signals  for  nitrogen  nuclei  in  /S-HMX  are  1.4  kHz  and  2.3 
kHz.  For  RDX  the  observed  linewidths  range  from  180  Hz  to  490  Hz.  The  room  temperature 
data  are  given  in  Table  10. 

Table  10.  '“‘N  NQR  data  for  nitrogen  atoms  of  jS-HMX  and  RDX. 

|8-HMX 


p_  =  3.737  MHz 

Ap  =  lA  kHz 

p_  =  3.623  MHz 

Ai'  =  2.3  kHz 

RDX 

p_  =  3.458  MHz 

Ap  -  0.49  kHz 

p_  =  3.410  MHz 

Ap  =  0.35  kHz 

p_  =  3.359  MHz 

Ap  =  0.43  kHz 

p+  =  5.24  MHz 

Ap  =  0.43  Hz 

p-\-  =  5.19  MHz 

Ap  =  0.18  Hz 

p+  =  5.04  MHz 

Ap  =  0.38  Hz 

The  contributions  to  the  second  moment  discussed  above  do  not  explain  the 
experimentally  observed  linewidths.  Thus  the  interactions  between  nitrogens  and  hydrogens  are 
considered.  In  the  absence  of  the  external  magnetic  field  there  is  no  contribution  to  the  second 
moment  from  the  N-H  dipolar  interaction  to  the  first  order.  However,  the  second  order 
contribution  may  be  significant  because  of  the  close  distances  between  nitrogen  and  hydrogen 
atoms  (Table  11  and  12). 

Table  11.  Data  for  /8-HMX.  The  distance  between  ring  nitrogen  atoms  and  the 
closest  hydrogen  atoms,  in  Angstroms. 


N(2)  - 

H(l) 

2.0854 

N(2)  - 

H(2) 

2.0870 

N(2)  - 

H(3’) 

2.1083 

N(2)  - 

H(4’) 

2.1145 

N(3)  - 

H(l) 

2.0671 

N(3)  - 

H(2) 

2.1181 

N(3)  - 

H(3) 

2.0937 

N(3)  - 

H(4) 

2.0479 
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Table  12.  Data  for  RDX.  The  distance  between  ring  nitrogen  atoms  and  the 
closest  hydrogen  atoms,  in  Angstroms. 


N(l)  -H(l) 

2.0781 

N(l)  -H(2) 

2.1038 

N(l)  -H(5) 

2.0849 

N(l)  -H(6) 

2.0977 

N(2)  -H(l) 

2.0372 

N(2)  -H(2) 

2.0869 

N(2)  -H(3) 

2.1036 

N(2)  -H(4) 

2.0643 

N(3)  -H(3) 

2.1027 

N(3)  -H(4) 

2.0598 

N(3)  -H(5) 

2.0475 

N(3)  -H(6) 

2.0841 

The  second  order  contributions  <w2>*'^(NH)  for  /S-HMX  are  presented  in  Table  13. 
Table  13.  Second  order  contributions  to  <w^>''^(NH)  for  ^-HMX. 


t 

N(2) 

N(3) 

rad/s 

rad/s 

0 

34 

35 

1 

175 

183 

2 

376 

392 

3 

624 

649 

However,  convergence  is  still  not  reached  for  t  —3.  The  computations  were  not  continued 
because  of  the  very  long  program  run  time.  The  values  of  <  are  much  larger  then 

<0)2  >1/2  calculated  in  the  first  order. 
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7.  SUMMARY 

This  work  focused  on  the  calculations  of  second  moments  of  '‘*N  NQR  signals  for  two 
explosive  materials,  jS-HMX  and  RDX.  Several  contributions  were  considered.  Computations 
included  interactions  between  completely  equivalent  nuclei,  interactions  between  nuclei  with  the 
EFG  tensors  of  the  same  magnitude  but  different  orientation,  and  the  second  order  contribution 
from  heteronuclear  dipolar  interactions. 

The  calculations  of  second  moments  to  the  first  order  show  that  <  for  /3-HMX 

is  larger  than  <  for  RDX  by  a  factor  of  two.  The  most  significant  difference  is  in  the 
contributions  from  the  completely  equivalent  nuclei,  which  are  almost  negligible  for  RDX. 

In  all  the  calculations  of  <  it  was  evident  that  the  distance  between  interacting 
nuclei  is  the  most  significant  factor.  The  angular  factor  influenced  the  value  of  SM  by  less  than 
one  order  of  magnitude  for  nuclei  separated  by  the  same  distance. 

The  numerical  results  of  the  calculations  of  second  moments  presented  in  this  work  do 
not  account  for  the  experimentally  observed  linewidths  of  NQR  signals.  However,  the 
second  order  contribution  from  heteronuclear  dipolar  interactions  for  j8-HMX  is  significant,  even 
though  the  convergence  was  not  reached.  Thus,  the  heteronuclear  dipolar  interactions  between 
nitrogens  and  hydrogens  may  be  important  factors  contributing  to  the  second  moments.  It  might 
be  confirmed  by  the  ‘‘‘N  NQR  study  of  deuterated  samples  of  the  explosive  materials. 
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ABBREVIATIONS 


a,  b,  c 

A 

E 

EFG 

g 

HMX 

/ 

NQR 

R 

RDX 


rf 

Q 

v» 


a,/3,7 

5, 


V 

y 

h 

(H 

X 

<Aa)2> 


unit  cell  parameters 
quadrupole  coupling  constant 
energy 

Electric  Field  Gradient 
metric  tensor 

cyclotetramethylene  tetranitramine 

angular  momentum 

Nuclear  Quadrupole  Resonance 

transformation  matrix 

cyclotrimethylene-trinitramine 

radio  frequency 

electric  quadrupole  moment 

principal  values  of  an  electronic  field  gradient  tensor 

unit  cell  parameters  (angles) 

Kronecker  delta 
Levi-Civita  symbol 
asymmetry  parameter 
magnetogyric  ratio 
Planck’s  constant 
Hamiltonian 
Lagrange  multiplier 
second  moment 
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APPENDIX  A 

The  transformation  matrix  between  the  coordinate  systems  defined  by  the  principal 
axes  of  EFG  tensors. 


In  this  appendix  we  derive  the  transformation  matrix  between  the  coordinate  systems 
defined  by  the  principal  axes  of  EFG  tensors  of  the  nuclei  in  the  magnetically  non-equivalent 
sites.  Let  vectors 


3 

j=l  > 


3 

1=1  » 


(41) 


denote  eigen-vectors  of  EFG  tensors  in  site  I  and  site  II.  Next  we  write  a  tensor 


E  = 


(42) 


where  the  elements  of  each  row  are  the  components  of  an  eigen-vector  of  the  EFG  tensor  in 
crystallographic  axes  system.  Similarly,  the  one-covariant  one-contravariant  tensor  Ffor  site  n 
is  constructed,  and  the  eigen-vectors  can  be  written  as 


3 


(43) 
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Further 


(Ej‘r'e,  =  (Ejr'Eta^, 

(£/)-■«,  = 

Oj  =  (£/)-«., 

/,  =  f/«y. 

/,  =  F/(£;)-‘e*. 

/,  =  (f  •£->),* 

«.  =  (£-f  ■‘)i,4 

R  =  EF 

Thus  the  matrix  R  transforms  the  base  vectors  e  into  /. 


(44) 
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APPENDIX  B 

TTie  transformation  of  x,  y,  z  coordinates  in  direct  lattice  to  the  orthonormal  system 
defined  by  the  principal  directions  of  an  EFG  tensor. 

First  we  define  a  vector  X  in  the  crystallographic  axes  system  c,- 

X  =  x^a^. 

Let  vector  Y  be  the  same  vector  but  defined  in  the  orthonormal  system  e,-  determined  by  the 
principal  directions  of  an  EFG  tensor. 

Y^y\. 

The  relation  between  two  sets  of  the  base  vectors  is 

(see  Appendix  A  for  the  definition  of  tensor  E).  If  vectors  X  and  Y  are  equal  then 

=  yU^y 

x^a^’a^  =  y^E^a^'a^ , 
x‘d{  =  y'Etsi, 
x>  =  y‘E{, 

(EjT  - 

=  (£/)-V. 

Thus  the  tensor  E  ‘^  transforms  coordinates  from  the  crystallographic  system  a,  to  the 
orthonormal  system  c,-. 
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APPENDIX  C 


Listing  of  programs 


%  ***  hmx_b.m  ***  Iolanta  M.  Snowel  January  26,  1995 
% 

%  This  program  contains  the  direct  unit  cell  axes,  fractional  atom  coordinates 
%  in  direct  lattice,  anisotropical  temperature  factors. 

%  It  calculates  the  reciprocal  metric  tensor  and  reciprocal  unit  cell  parameters. 

%  In  general  calculations  are  performed  in  a  non-orthogonal  basis. 

% 

%  g  -  direct  metric  tensor 
%  G  -  reciprocal  metric  tensor 

%  a,  b,  c,  alpha,  beta,  gamma  -  direct  unit  cell  parameters 

%  ar,  br,  cr.  Alpha,  Beta,  Gamma  -  reciprocal  unit  cell  parameters 

% 

%  crystal:  HMX_B 
% 

%  Chang  S.  Choi,  Henry  P.  Boutin 
%  Acta  Cryst.  1970,  B26,  1235-1240. 

a=6.54;  b= 11.05;  c=8.7;  alpha=pi/2;  beta = (124. 3/360)*(2*pi);  gamma=pi/2; 
%beta=(124.3/360)*(2*pi); 
g=[  a*a  a*b*cos(gamma)  a*c*cos(beta) 
b*a*cos(gamma)  b*b  b*c’^cos(alpha) 
c*a*cos(beta)  c*b*cos(alpha)  c*c  ]; 


G=inv(g); 

ar=sqrt(G(l,l));  br=sqrt(G(2,2));  cr=sqrt(G(3,3)); 


Alpha=acos(G(2,3)/(sqrt(G(2,2)*sqrt(G(3,3))))); 
Beta=acos(G(l  ,3)/(sqrt(G(l ,  l)*sqrt(G(3,3))))); 
Gamma =acos(G(  1 ,2)/(sqrt(G(  1 ,  l)*sqrt(G(2 ,2))))) ; 


data  =  [ 


-5928 

-21 

-2920 

1092  504 

651  61 

148  242 

-3482 

-249 

-2058 

1199  363 

617  -62 

281  7 

-188 

1228 

-377 

1366  277 

566  -26 

482  -75 

-1092 

2023 

316 

1717  279 

834  52 

815  -39 

-6885 

735 

-4182 

1602  600 

1049  190 

73  2 

-6991 

-598 

-2372 

1708  823 

1282  -237  1121  -104 

422 

2460 

1818 

3062  585 

1361  -53 

1155  -456 

-3290 

2206 

-609 

1883  603 

1200  398 

914  28 
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-1967 

660 

-2158 

2448  1 

147 

547 

-879 

218 

-2623 

-3146 

1326 

-3223 

3322 

2031 

1130 

2793 

867 

-494 

]; 


1103  357  467  -27  456  -57 
1865  247  755  -102  840  -31 
3247  585  1081  27  1623  -183 
2320  818  906  -231  383  207 
1478  636  1961  80  310  -111 
3518  1199  1008  212  1520-1 


xyz=data(:,l:3);  xyz =0.0001  *xyz; 


%  B  -  thermal  parameter  tensor.  It  must  be  transformed  to  beta  form: 

% 

%  exp(-(Bll*h^2  +  Bll*Vr2  +  B33*r2  +  2*B12*h*k  +  2*B13*h*l  +2*B23*k*l)) 


Bm  =0.00001  *data(:,  4:9);  clear  data; 


Bm(:,4:6)  =  Bm(:,4:6)/2; 


i|e:((3|e:|es|c3ie3|e3|e3|e3|c:|c3Ka|(  3<c3|c3te]|(9|e:(e3|e « 3|e:(c^  a|c  #3(e3(C3|(9(e:(c9|e  s(e3|e]|c3|cs|c:|o(c3|e9KaK3ie:|e  ife** 

3|e  3te  3|c  ]|e  4e  ^  3(c  3|e  Ke  3|e  He  :|c  #  a|c  3|c  a(c  :|e « Die  ♦  9|e  9K  a|c  9|e  3((  3(e  « 9|e  3K  3|c  3(c  #  3|e  )|c  a|o|e  :(e  9K  afe  sK  >|e  aK  3(e  aK  3K  3tc 

%  ♦**  rdx.m  **♦  Jolanta  M.  Snowel  May  15,  1995 
% 

%  This  program  contains  the  direct  unit  cell  axes,  fractional  atom  coordinates 
%  in  direct  lattice,  anisotropical  temperature  factors. 

%  It  calculates  the  reciproc^  metric  tensor  and  reciprocal  unit  cell  parameters. 
%  In  general  calculations  are  performed  in  a  non-orthogonal  basis. 

% 

%  g  -  direct  metric  tensor 
%  G  -  reciprocal  metric  tensor 

%  a,  b,  c,  alpha,  beta,  gamma  -  direct  unit  cell  parameters 

%  ar,  br,  cr.  Alpha,  Beta,  Gamma  -  reciprocal  unit  cell  parameters 

% 

%  crystal:  RDX 
% 

%  Chang  S.  Choi,  E.  Prince 
%  Acta  Cryst.  1972,  B28,  2857-2862. 

a=13.182;  b=11.574;  c=10.709;  alpha=pi/2;  beta=pi/2;  gamma=pi/2; 

Z_au  =  8;  %  number  of  asymmetric  units 

g=[  a*a  a*b*cos(gamma)  a^c^cosCbeta) 
b*a*cos(gamma)  b*b  b*c*cos(alpha) 
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c’^a’'‘cos(beta)  c*b*cos(alpha)  c*c  ]; 


G=inv(g); 

ar=sqrt(G(l,l));  br=sqrt(G(2,2));  cr=sqrt(G(3,3)); 


Alpha=acos(G(2,3)/(sqrt(G(2,2)*sqrt(G(3,3))))); 
Beta = acos(G(  1 , 3)/(sqrt(G(  1,1)  *sqrt(G(3 , 3))))) ; 
Gamma=acos(G(l,2)/(sqrt(G(l,l)*sqrt(G(2,2))))); 

data  =  [ 

1839  3578  4400  322  297  206  -70  -69  7 

503  2440  3395  233  145  344  8  -1  -17 

1487  3813  2159  252  332  175  -34  34  -31 

1761  4360  3330  267  266  251  -109  1  0 

877  2996  4537  295  210  250  14  33  9 

536  3222  2329  211  228  274  15  -35  -37 

2260  5379  3346  217  241  348  -55  79  -26 

155  3525  5294  455  322  256  157  104  77 

-333  3875  2078  284  306  316  72  -72  -53 

2270  5932  2388  566  330  407  -122  82  68 

2649  5687  4338  342  330  417  -86  -13  -131 

-693  3142  5262  333  718  426  138  153  156 

454  4271  5990  894  381  258  251  57  -53 
-1121  3534  2504  246  572  584  44  -62  22 

-236  4691  1390  644  441  415  185  -91  62 

2401  2944  4226  372  532  564  -14  -128  141 

2013  4068  5244  762  560  270  -201  -127  -59 

-261  2117  3540  371  433  597  -146  -18  16 

1015  1728  3192  464  263  594  89  -16  -114 

2052  3166  1917  314  602  524  46  80  -255 

1440  4432  1413  637  602  363  -170  17  108 

]; 


xyz=data(:,l:3);  xyz=0.0001*xyz; 

%  B  -  thermal  parameter  tensor.  It  must  be  transformed  to  beta  form: 

% 

%  exp(-(Bll*h"2  +  B22*k^2  +  B33*r2  -I-  2*B12*h*k  -I-  2*B13*h*l  -F2*B23*k*l)) 
Bm  =0.01*data(:,4:9)/4; 


Bm(:,l)  =  Bm(:,l)*G(l,l); 
Bm(:,2)  =  Bm(:,2)*G(2,2); 
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Bin(:,3)  =  Bm(:,3)*G(3,3); 
Bm(:,4)  =  Bm(:,4)*ar*br; 
Bni(:,5)  =  Bm(:,5)*ar*cr; 
Bm(:,6)  =  Bm(:,6)*br*cr; 


atom_name  =  [’C(1)’;’C(2)’;’C(3)’; 

’N(l)’ ;  ’N(2)  ’ ;  ’N(3)’ ;  ’N(4)’ ;  ’N(5)’ ;  ’N(6)’ ; 
‘0(1)’;’0(2)’;’0(3)’;’0(4)’;’0(5)’;’0(6)’; 
’H(1)’;’H(2)’;’H(3)’;’H(4)’;’H(5)’;’H(6)’  ]  ; 

♦♦♦♦*♦#♦*♦♦♦»#♦♦#♦*♦♦♦#**»*♦♦♦*♦»*♦*♦♦♦*»♦♦♦♦♦***♦**»***»♦****♦**** 


%  ♦**  bonds. m  ***  Jolanta  M.  Snowel  January  27,  1995. 

%  This  program  calculates  distances  between  atoms. 

%  Calculations  are  performed  in  a  non-orthogonal  basis. 

%  g  -  the  direct  metric  tensor 
%  G  -  the  reciprocal  metric  tensor 
%  a,  b,  c,  alpha,  beta,  gamma  -  unit  cell  parameters 
%  ar,  br,  cr  Alpha,  Beta,  Gamma  -  reciprocal  unit  cell  parameters 

%atoms  %  ♦.m  file.  It  loads  data:  unit  cell  parameters,  x,  y,  z,  and  Bij  (beta). 
hmx_b 

%  g=[  a*a  a*b*cos(gamma)  a*c*cos(beta) 

%  b^a^cosCgamma)  b*b  b*c*cos(alpha) 

%  c’'‘a*cos(beta)  c*b*cos(alpha)  c*c  ] 


%  G  =  inv(g); 


%  BOND  LENGTH  (or  distance  between  2  atoms) 

k  =  size(xyz);  k  =  k(l);  1=0; 

for  i  =  l:k 

for  j  =  (i+l):k 

V  =  xyz(i,:)-xyz(j,:); 

d  =  sqrt(v*g*v’); 

%  ifd  <  1.5 

1=1+ 1; 

distance(l,:)  =  [  i  j  d  ]  ; 

%  end 

end 
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end 


fprintf(’\n  \n  Bond  lengths  in  Angstroms  \n  \n’,  distance’); 
fprintf(’%6.0f  -  %3.0f  %  12. 4f\n’, distance’); 

^  \q  3^  ♦♦*3*c**3jc* J|«*3fcj|c3(c:}c>|ej»ca|e>|c 

fid  =  fopen(’c:\matlab\jolanta\crystal\bond.out’,’w’); 

fprintf(fid,’DATA  FOR  beta-HMX.\n  \n  ’); 

fprintf(fid,’\n  \n  Fractional  coordinates  in  the  direct  lattice  \n  \n’); 

fprintf(fid,’%12.4f  %12.4f  %12.4f\n’,  xyz’); 

fi)rintf(fid,’\n  \n  Bond  lengths  in  Angstroms  \n  \n’); 

f^rintf(fid,’%6.0f  -  %3.0f  \t  %3.4An’,  distance’); 

fclose(fid) 

3|c  3ie  sic  3^  sf;  3^  9|e  9K  3k  :|£ 3)c  lie  ♦  aK  *  3f:  :(c  a|c  *  afe  af:  3f^  9^ 

Ik  ak  ik  afe  ak  ak  3|c  ak  ife  a(e  afe  ak  3|e  a(e  afe  ak  3k  ak  ak  Ik  3{e  ak  3k  3k  ak  ik  ak  ak  ak  ak  ak  ak  ik  ak  :k  ak  ak  ak  ak  ak  ak  a|e  ate  ak  :ie  )k  Ik  ak  ak  ak  ak  ak  ak  ak  ak  ak  ak 


%  *♦*  angles. m  ♦**  January  27,  1995. 

%  This  program  calculates  angles  between  bonds. 

%  It  requires  input  file  such  as  HMX_B.M 

%  Make  a  table  of  atom  numbers: 


5  1  6 

2  1  5 

2  1  6 

1  2  9 

2  9  3 

2  9  11 

2  9  12 

3  9  11 

3  9  12 

11  9  12 

7  4  8 

3  4  7 

3  4  8 

4  3  9 

4  3  10 

9  3  10 

3  10  14 

3  10  13 

13  10  14  ]; 
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%  ANGLES 


k  =  size(at);  k  =  k(l);  ang  =  zeros(k,l); 


for  i  =  1  :  k 

vecl  =  [  at(i,2)  at(i,l)]; 

vecl  =  xyz(vecl(l),:)  -  xyz(vecl(2),:); 

vec2  =  [  at(i,2)  at(i,3)  ]; 

vec2  =  xyz(vec2(l),:)  -  xyz(vec2(2),:); 

vecll=  sqrt(vecl*g*vecr); 

vec21=  sqrt(vec2*g*vec2’); 

ang(i)  =  57.2958*acos(vecl*g*vec2V(vecll*vec21)); 
end 

bond_angles  =  [at  ang]; 

^rintf(’\n  \n  Bond  angles  for  beta-HMX.\n  \n’); 
Q)rintf(’%12.0f  -%3.0f  -%3.0f  \t  %12.1An’,  bond_angles’); 


fid  =  fopen(’c:\matlab\jolanta\crystal\angles.out’,’w’); 
fprintf(fid,’DATA  FOR  beta-HMX.\n  \n  ’); 

^rintf(fid,’\n  \n  Bond  angles  for  beta-HMX.\n  \n’); 

^rintf(fid,’%12.0f  -%3.0f  -%3.0f  \t  %12.1f\n’,  bond_angles’); 
fclose(fid) 

ititu***************************************************************** 


%  quadric.m  January  27,  1995. 

%  This  program  calculates  principal  axes  of  thermal  ellipsoids 
%  and  eigenvectors.  The  eigenvectors  are  orthogonal  in  direct  basis. 
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%  B  is  a  matrix  of  beta  parameters 


clear  all; 
hmx_b 

m  =  size(Bm);  m  =  m(l);  j  =  0;  d  =  zeros(m,3);  v  =  zeros(3,3*m); 
for  k  =  l:m 

B(l,l)  =  Bm(k,l);  B(2,2)  =  Bm(k,2);  B(3,3)  =  Bm(k,3); 

%B(1,2)  =  Bm(k,4)/2;  B(l,3)  =  Bm(k,5)/2;  B(2,3)  =  Bm(k,6)/2; 
B(l,2)  =  Bm(k,4);  B(l,3)  =  Bm(k,5);  B(2,3)  =  Bm(k,6); 

B(2,l)  =  B(l,2);  B(3,l)  =  B(l,3);  B(3,2)  =  B(2,3); 


Bg=(B*g); 

[V,D]=eig(Bg); 

d_aux = (((diag(D))/(2*pi’^pi)) .  "0.5)  * 1000; 

%d=(((diag(D))/(2*pi*pi))."0.5); 

d(k,:)=d_aux’; 

v(l:3J+l:j+3)=V’;  j  =  j  +  3; 
end 

d  =  d’  ;  d  =  sort(d);  d  =  d’;  clear  d_aux; 

fprintf(’\n  \n  Magnitudes  of  the  principal  axes  of  the  thermal  ellipsoids.  \n  \n’); 
fprintfC%12.1f  %12.1f  %12.1f  \n’,  d’); 

fiprintf(’\n  \n  Eigenvectors  of  the  thermal  ellipsoids.  \n  \n’); 

fprintfC %  12.4f %  12.4f%  12.4f\n%  12.4f%  12.4f%  12.4f\n%  12.4f%  12.4f%  12.4An\n’,  v); 
^  *^f^f********  Writing  eigenvalues  and  eigenvectors  to  QUADRIC. OUT 

]|i «  *  « i|c He 

fid  =  fopen(’c:\matlab\jolanta\crystal\quadric.out’,’w’); 

^rintf(fid,’  DATA  FOR  beta-HMX.  \t\t\t\t\t\t\t  January  27,  1995.’); 

^rintf(fid,’\n  \n  \n  \n  Thermal  parameters  \n  \n’); 

fprintf(fid,’%4.0At%4.0At%4.0At%4.0At%4.0f\t%4.0At\n’,  (data(:,4:9))’  ); 
fprintf(fid,’\n  \n  Magnitudes  of  the  principal  axes  of  the  thermal  ellipsoids.  \n  \n’); 
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fprintf(fid,’%12.1f  %12.1f  %12.1f  \n’,  d’); 

Q)riiitf(fid,’\An  \n  Eigenvectors  of  the  thermal  ellipsoids.  \n  \n’); 

fprintf(fid,’  %  12.4f%  12.4f%  12.4f\n%  12.4f%  12.4f%  12.4An%  12.4f%  12.4f%  12.4f\n\n’, 

V); 

fclose(fid) 

♦  ♦♦♦♦UtKc*#**************************’''******************************* 


%  ♦**  sec_mom.m  ***  Jolanta  M.  Snowel  June  29,  1995, 

% 

%  This  program  calculates  contributions  to  the  second  moment  <w‘^2> 

%  from  COMPLETELY  EQUIVALENT  NUCLEI  and 
%  from  EQUIVALENT  NUCLEI  with  Different  Directions  of  EFG  Tensors. 

%  Nitrogen  nuclei:  spin  1=1,  asymmetry  parameter  !=  zero, 

%  no  magnetic  freld. 

%  LAST  Version  :  sec_mom.j21,  June  21,  1995. 

%  Two  sites  are  considered  : 

%  I  -  N(2)  in  the  positions:  (x,y,z)  &  (-x,-y,-z) 

%  n  -  N(2)  in  the  positions:  (-x,y-l/2,-z-l/2)  &  (x,-y+l/2,z+l/2) 

%  (-x,y-l/2,z-l/2)  =  -(x,-y+l/2,z+l/2) 

%  <w'"2>  eqs  :  2,15  2.16  &  2.33  S.Vega 

clear  all 

czas  =  cputime; 

hmx_b  %  hmx_b.m  file.  It  loads  data:  unit  cell  parameters,  x,  y,  z,  and  B_ij  (beta). 

clear  Alpha  Beta  Gamma  alpha  beta  gamma  ar  br  cr  a  b  c  Bm; 

t  =  3;  %  Translation  limits  for  n,  m,  and  p. 

atom  =  2;  %  1  for  N(2)  ;  2  for  N(3) 

=  fopen(’c:\matlab\jolanta\crystal\hmx_tab.out’,’a’); 

%  for  atom  =  [12] 

fprintf(fp,’Data  for  HMX;  \t  N(%1.0f):\t  translation  limits  t  =  %1.0f;  \n  \n’,atom+l,t); 
^rintf(fp,’\t  cosi  cosll  <w2>a,II)  <w2>a,I)  <w2>(n,n)  <w2>-l/2 

\n  \n’); 

Q)rintf(’Data  for  HMX;  \t  N(%1.0f):\t  translation  limits  t  =  %1.0f;  \n  \n’,atom+l,t); 
:fprintf(’\t  cosi  cosH  <  w2  >  (I,n)  <  w2  >  (1,1)  <  w2  >  (n,n)  <  w2  >  -1/2  \n  \n’); 
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^  j|ca|ej|c3*e**^*#J|c*JKJ|es|c*3|ej|e3jcj(e^JlcHe  gpQ  AXES 


%  Make  a  table  of  pairs  of  atom  numbers  in  each  row, 

%  numbers  refer  to  the  sequence  in  ’xyz’  matrix  in  hmx_b.m  file 
%  (also  see  efg.m,  efg_uc.m,  &  efg_ucat.m): 


at  =  [ 

6  5  12  %  EFG  for  N(2) 

8  7  4  3  ];  %  EFG  for  N(3) 

at  =  at(atom,;); 

sqrt_g  =  sqrt(det(g)); 

k  =  size(at,l);  efg  =  zeros(k*12,3);  RR  =  zeros(k*12,3); 

vector!  =  zeros(4*k,3);  vector2  =  zeros(4*k,3);  vector3  =  zeros(4*k,3); 


for  i  =  1  ;  k 

vecl  =  xyz(at(i,2),:)  -  xyz(at(i,l),:);  %  [EFG]1  =  [O  -  O]  < —  auxiliary 

vec2  =  xyz(at(i,4),:)  -  xyz(at(i,3),:);  %  [EFG]2  =  [N  -  N] 

vecll=  sqrt(vecl*g*vecr); 
vec21=  sqrt(vec2*g*vec2’); 

vecl  =  vecl/vecll; 
vec2  =  vec2/vec21; 

vector2(i,:)  =  vec2; 

^  »**]K*)ic*««**  p£pG]3  =  [O  -  O]  X  [N  -  N]  ****’t‘******» 

vec3(l)  =  sqrt_^*(  vecl(2)*vec2(3)  -  vecl(3)*vec2(2) )  ; 
vec3(2)  =  sqrt_^*(  vecl(3)*vec2(l)  -  vecl(l)*vec2(3) )  ; 
vec3(3)  =  sqrt_g*(  vecl(l)*vec2(2)  -  vecl(2)*vec2(l) )  ; 

vector31  =  sqrt(vec3*G*vec3’); 
vector3(i,:)  =  (vec3/vector31)*G  ; 
vec3  =  vec3*G  ; 

4,**********  [EFGjl  =  [EFG]2  x  [EFG]3  ************* 
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vectorl(i,l)  =  sqrt _g*(  vec2(2)*vec3(3)  -  vec2(3)*vec3(2) )  ; 
vectorl(i,2)  =  sqrt_g*(  vec2(3)*vec3(l)  -  vec2(l)*vec3(3) )  ; 
vectorl(i,3)  =  sqrt^*(  vec2(l)*vec3(2)  -  vec2(2)*vec3(l) ) ; 

vectorl(i,:)  =  vectorl(i,:)*G  ; 

vectorll  =  sqrt(vectorl(i,:)*g*vectorl(i,:)’); 

vectorl(i,:)  =  vectorl(i,:)/vectorll; 

end  %  *♦*♦*♦**♦**♦  END  *♦*♦♦*****♦**♦ 

%  eFG  axes  of  all  N(2)  and  N(3)  in  equivalent  positions  ******* 

vectorl(k+l:2*k,:)  =  [-vectorl(l:k,l)  vectorl(l:k,2) -vectorl(l:k,3)];  %  2) 

-x,y+l/2,-z+l/2 

vector2(k+l:2*k,:)  =  [-vector2(l;k,l)  vector2(l:k,2)  -vector2(l:k,3)];%  +  0.5 
vector3(k+l:2*k,:)  =  [-vector3(l:k,l)  vector3(l:k,2)  -vector3(l:k,3)];%  +  0.5 


vectorl(2*k+l:3*k,:)  = -vectorl(l:k,:)  ;  %  (3) -x,-y,-z 

vector2(2*k+l:3*k,:)  =  -vector2(l:k,:)  ; 

vector3(2*k+l:3*k,:)  =  -vector3(l;k,:)  ; 

vectorl(3*k+l:4*k,:)  = -vectorl(k+l:2*k,:);  %  (4)  x,-y+l/2,z+l/2 

vector2(3*k+l:4*k,:)  =  -vector2(k+l:2*k,:);  %  +  0.5  ; 

vector3(3*k+l:4*k,:)  =  -vector3(k+l:2*k,:);  %  +  0.5  ; 

^  3iei|ej|e«a(e:(C3|e3|ea|ea|e9(c 

E  =  [vectorl(l,:)  ;  vector2(l,:)  ;  vector3(l,:)];  invE  =  inv(E)  ; 

%  E  =  [vectorl(3,:)  ;  vector2(3,:)  ;  vector3(3,:)];  invE  =  inv(E)  ; 

n=l;  1=3; 

for  i  =  1  :  4*k 

efg(n:l,:)  =  [vectorl(i,:)  ;  vector2(i,:)  ;  vector3(i,:)]; 

F  =  efg(n:l,:); 

RR(n:l,:)  =  E*inv(F); 
n  =  1+  1;  1  =  1  +  3; 
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end 


%  »♦♦♦*♦♦♦*♦♦  end  of  calculations  of  the  efg  tensors 

#  ate  3|e  ^  ate  *  ♦  3(c  :tc  :(c 

%  ***********  &  TRANSFORMATION  MATRIX  [RR] 

ate  ate  a|e  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  at:  at(  ate  ate  3tc  ate  at:  ate  ate  ate  ate 


% 

ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  a|c  ate  ate  ate  ate  at:  at:  ate  ate  ate  ate  ate  ate  ate  aK  ate  ate  ate  a|e  ate  at:  ate  ale  ate  ate  ate  ate  a|c  ate  ate  ate  ate  ate  at:  ate  ale  ate  ate  ale  ate  ate  ate  ate  a|e  ate  a^ 

ate  ate  ate  ate  ate  ate  ate 

%  ♦♦**♦*♦♦♦♦  POSITIONS  OF  EQUIVALENT  NUCLEI  {  N(2)  } 

ate  ate  a|e  ate  ate  ate  ate  ate  ate  ate  at:  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate  ate 

%  Either  atoms  N(2)  or  N(3)  from  ’xyz’  matrix  (hmx_b.m  file)  are  considered: 
m  =  1  ; 

xyz_d  =  zeros(m*4,3); 

xyz_d(l,:)  =  xyz(atom+l,:);  %  fractional  coordinates  of  N(2)  or  N(3) 


%  *♦*  xyz_d  -  matrix  of  the  coordinates  of  nitrogen  N(2)  atoms 
%  ***  in  the  following  positions: 

%  ***  (1)  x,y,z  (2) -x,y-l/2,-z-l/2  (3) -x,-y,-z  (4)  x,-y+l/2,z+l/2  . 

%  ***  (1)  &  (3)  -  1st  molecule;  (2)  &  (4)  -  2nd  &  3rd  molecule; 

xyz_d(m+l:2*m,l)  =  -xyz_d(l:m,l)  ;  %  2)  -x,y-l/2,-z-l/2 

xyz_d(m+l:2’'‘m,2)  =  xyz_d(l:m,2)  +  0,5  ; 

xyz_d(m+l:2*m,3)  =  -xyz_d(l:m,3)  +  0.5  ; 


%xyz_d(m+l:2*m,l)  = -xyz_d(l:m,l)  ;  %  2) -x,y+l/2,-z+l/2 

%xyz_d(m+l:2*m,2)  =  xyz_d(l:m,2)  -  0.5  ; 

%xyz_d(m+l:2*m,3)  =  -xyz_d(l:m,3)  -  0.5  ; 

xyz_d(2*m+l:3*m,l)  = -xyz_d(l:m,l)  ;  %  (3) -x,-y,-z 

xyz_d(2’'‘m+l:3*m,2)  =  -xyz_d(l:m,2)  ; 

xyz_d(2*m+l:3*m,3)  =  -xyz_d(l:m,3)  ; 


xyz_d(3*m+l:4*m,l)  =  xyz_d(l:m,l)  ;  %  (4)  x,-y+l/2,z+l/2 

xyz_d(3*m+l:4*m,2)  =  -xyz_d(l:m,2)  +  0,5  ; 

xyz_d(3*m+l:4*m,3)  =  xyz_d(l:m,3)  +  0.5  ; 
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***«»*«  (^sines  ***♦♦•***** 


k  =  size(efg,l);  countl  =  0;  count2  =  0; 
i=l; 

forj  =  l:(k/3) 
vic2(j,:)  =  efg(i,:)*invE; 
i=i+3; 
end 

^  4:3|es|ea|c3|e]|c#3|ei(e:|e4e3|ea|csie  3|e:|c 

%Omega  =  zeros(216,4);  ANGl  =  zeros(216,l);  ANG2  =  zeros(216,l);  DA  = 
zeros(216,l); 

%CosI  =zeros(216,l);  CosH  =zeros(216,l); 

count2  =  0;  count3  =  0; 

%cl=6;  c2=2;  c3=3;  [-5  -3  -1  1  4  11]  ;  [-3  -1  5] 

%for  cl  =  [-11-24] 

%for  c2  =  [-4  -1  3] 

%for  c3  =  [  -3  0  5] 

for  cl  =4 
for  c2  =3 
for  c3  =0 

vicl  =  [cl  c2  c3]/(sqrt([cl  c2  c3]*[cl  c2  c3]’)); 
for  i  =  1:4 

cos_table(i)  =  [vicl*vic2(i,:)’]; 
end 

countl  =countl  -I- 1 ; 
cosinus_table(countl,:)  =  cos_table; 

^  »*««**itc«**  SECOND  MOMENT  CONTRIBUTION  FROM  EQUIVALENT  NUCLEI 
^  DIFFERENT  DIRECTIONS  OF  THE  EFG  TENSORS 

3|e3|e:|e#9|e:|e3|e4e4e4e« 

SecMoment  =  zeros(6,l);  SecMom  =  0;  %SecondMom  =  0; 

Da  =  0;  Angl  =  0;  Ang2  =  0; 
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^  DistEHCC  bCtWCCIl  2  EtOITlS  ♦******’*«****’*^******** 


k  =  size(xyz_d,l);  1=0; 

for  i  =  1:4 

cosi  =  cos_table(i); 
vl  =  xyz_d(i,:); 
for  h  =  0: 1 

h  =  i  +  2*h; 
j  =  rem(h,4)  +1; 

cosll =cos_table(j); 

R  =  RR(3*j-2:3*j,:);  %  ***  Transformation  matrix  **** 

v2  =  xyz_dO,:); 
for  m  =  -t:t 

v2(l)  =  xyz_d(j,l)  +  m  ; 

for  n  =  -t:t 

v2(2)  =  xyz_d(j,2)  +  n  ; 
for  p  =  -t:t 

v2(3)  =  xyz_d(j,3)  +  p  ; 

V  =  v2  -  vl  ; 
r  =  sqrt(v*g*v’)  ; 

1=1+1; 

D  =  v*invE  ; 

%  d_o  =  sqrt(I)*D’)  ; 


X  =  D(l)/r;  Y  =  D(2)/r;  Z  =  D(3)/r; 

XOS  =  1  -  3»D(l)»D(l)/(r»r); 

YOS  =  1  -  3*D(2)*D(2)/(r»r); 

ZOS  =  1  -  3*D(3)*D(3)/(r*r); 


count2  =  count2  +  1; 

%  ***  Dipolar  Hamiltonian  *** 

%  hi  =  h_23  =  h_32;  h2  =  h_56  =  h_65;  h3  =  h_78  =  h_87; 

hi  =  R(3,3)*ZOS  -  3*R(2,3)*Y*Z  -  3*R(1,3)*X*Z  ; 
h2  =  R(l,l)*XOS  -  3»R(2,1)*X*Y  -  3*R(3,1)*X*Z  ; 


44 


h3  =  R(2,2)*YOS  -  3*R(1,2)*X*Y  -  3*R(3,2)*Y*Z  ; 


%  *♦*  SECOND  MOMENT  for  Equivalent  Nuclei  with  Different  Directions  of  the  EFG 
Tensors. 

%  ***  eq.  2.33  ♦** 

angl  =  cosI*cosn/(cosI*cosI  +  cosII*cosII); 

dipol  =  2*hl*hl  +  2*h3*h3  +  4*h2*h2  -  8*hl*h3*angl; 

SM  =  dipol/(r*r*r*r*r*r); 

%SM  =  ( 2*hl*hl  +  2*h3*h3  +  4*h2*h2  -  8*hl*h3*cosI*cosn/(cosI*cosI  + 
cosn*cosII))/r‘^6; 

SecMom  =  SecMom  +  SM  ;  %  Da  =  Da  +  da;  Angl  =  Angl  + 

angl; 

%  Distanced, :)  =  [ijninpr  SM  ]  ; 

%  angl  =  0; 

end  %  end  of  ’  for  p  =  -t:t  ’ 
end  %  end  of  ’  for  n  =  -t:t  ’ 
end  %  end  of  ’  for  m  =  -t:t  ’ 
end  %  end  of  ’  for  h  =  0: 1  ’ 

SecMomentfi)  =  SecMom; 

SecMom  =  0  ; 

SM  =  0; 

end  %  end  of  ’  for  i  =  1:4  ’ 

%  DA(countl)  =  Da/8;  Da  =  0; 

%  ANGl(countl)  =  Angl/8; 

SecMoment  =  SecMoment*2592200;  %  const  =  (l/6)*gamma'^4*h_bar''2 

%fprintf(’\n  \n  Distance  between  2  atoms  in  Angstroms  and  Second  Moments  \n  \n’); 
%fprintf(’\n  \n  Equivalent  nuclei  with  different  directions  of  EFG  tensors\n  \n’); 
%4rintf(’%6.0f  -  %3.0At  %3.0f  %3.0f  %3.0f  %12.4f%24.8f  \n’, Distance’); 

%^rintf(fp,’\n  \n  Distance  between  2  atoms  in  Angstroms  and  Second  Moments  \n  \n’); 
%^rintf(Q),’\n  \n  Equivalent  nuclei  with  different  directions  of  EFG  tensors\n  \n’); 
%^rintf(^,’%6.0f  -  %3.0f\t  %3.0f  %3.0f  %3.0f  %12.4f%24.8f  \n’, Distance’); 
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%  ♦♦***♦♦♦**  SECOND  MOMENT  CONTRIBUTION  FROM  COMPLETELY  EQUIVALENT 
NUCLEI  *********** 

SecMom  =  0; 

%  **********  Distance  between  2  atoms  ********************* 

k  =  size(xyz_d,l);  1=0;  %  Distance  =  zeros((k*k-k)/2,4); 

cosi  =  cos_table(l);  cosll  =cos_table(2); 

for  i  =  1:2 

vl  =  xyz_d(i,:); 
forj  =  i:2:i+2 

vie  =  xyz_d(j,:); 


for  m  =  -t:t 

vle(l)  =  xyz_d(j,l)  +  m; 
for  n  =  -t:t 

vle(2)  =  xyz_d(j,2)  +  n; 
for  p  =  -t:t 

vle(3)  =  xyz_d(j,3)  +  p; 

V  =  vie  -  vl  ; 
r  =  sqrt(v*g*v’)  ; 
if  r  >  0.1 
1=1+1; 

D  =  v*invE  ; 

%  d_o  =  sqrt(D*D’)  ; 

X  =  D(l)/r;  Y  =  D(2)/r;  Z  =  D(3)/r; 

XOS  =  1  -  3*D(l)*D(l)/(r*r); 

YOS  =  1  -  3*D(2)*D(2)/(r*r); 

ZOS  =  1  -  3*D(3)*D(3)/(r*r); 
count3  =  count3  +  1; 

%  <w"2>  =  sum(  8*XOS^2  +  4’^YOS"2  +  4*ZOS"2  -  8*YOS*ZOS)/r"6  % 
completely  equivalent  nuclei 

%  *♦*  SECOND  MOMENT  for  Completely  Equivalent  Nuclei.  ******* 

%  ***  eq.  2.15  **♦ 
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SM  =  (8*X0S*X0S  +  4*Y0S*Y0S  +  4*Z0S*Z0S  -  8*Y0S*Z0S)/(r*r*r*r*r*r); 


SecMom  =  SecMom  +  SM  ; 

%  *♦*  Dipolar  Hamiltonian  ♦** 

%  Distance(l,:)  =  [  i  j  m  n  p  r  SM  ]  ; 

end  %  end  of  ’  if  r>0.1  ’ 
end  %  end  of  ’  for  p  =  -t:t  ’ 

end  %  end  of  ’  for  n  =  -t:t  ’ 

end  %  end  of  ’  for  n  =  -t:t  ’ 

end  %  end  of  ’  for  j  =  i:2:i+2  ’ 

ang2  =  cosI*cosI/(cosI*cosI  +  cosn*cosII); 

SecMoment(i+4)  =  1296100*SecMom*ang2;  %  const  = 

(l/12)*gammaM*h_bar^ 

SecMom  =  0  ; 

%fprintf(’  \n  Distance  between  2  atoms  in  Angstroms  and  Second  Moments  \n  \n’); 
%fprintfC  \n  Equivalent  nucleiXn  \n’); 

%fprintf(’%6.0f  -  %3.0f\t  %3.0f  %3.0f  %3.0f  %12.4f%24.8f  \n’, Distance’); 

%fprintf(fp,’  Distance  between  2  atoms  in  Angstroms  and  Second  Moments  \n  \n’); 
%fprintf(fp,’\n  Equivalent  nuclei\n  \n’); 

%^rintf(fp,’%6.0f  -  %3.0At  %3.0f  %3.0f  %3.0f  %12.4f%24.8f  \n’, Distance’); 
end  %  end  of  ’  for  i  =  1:2  ’ 

Omega_2(l)  =  (sum(SecMoment(l:4)))/4; 

Omega_2(2)  =  SecMoment(5)/2; 

Omega_2(3)  =  SecMoment(6)/2; 

Omega_2(4)  =  sqrt((sum(SecMoment(l:4)))/4+0.5*SecMoment(5)+0.5*SecMoment(6)); 

Omega(countl,:)  =  Omega_2; 

SecondMom(countl,:)  =  SecMoment’; 

omega_2(l)  =  (sum(SecondMom(l:4)))/4; 
omega_2(2)  =  SecondMom(5)/2; 
omega_2(3)  =  SecondMom(6)/2; 
omega_2(4)  = 

sqrt((sum(SecondMom(l:4)))/4+0.5*SecondMom(5)+0.5*SecondMom(6)); 

omega_2m  =  mean(Omega); 
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fprintf(^,’\t  %  10.3f  %  10.3f  %  12.0f%  12.0f%  12.0f  %  14.0f 

\n’  ,cosI,cosn,Oniega(countl , 

fprintf(’\t  %  10.3f  %  10.3f  %  12.0f%  12.0f %  12.0f  %  14.0f 

\n’,cosI,cosn,Omega(countl , 

end  %  end  of  ’  for  c3  =  ’ 

end  %  end  of  ’  for  c2  =  ’ 

end  %  end  of  ’  for  cl  =  ’ 

^rintf(lp,’\n\n’); 

f^rintf(’\n\n’); 

%fprintf(’ %8.3At%8.3At%  10.0At%  10.0At%  10.0At%  lO.Of  \n’,[cosinus_table(: ,1:2) 
Omega]’); 

%fprintf(Q),’\n\n%8.3At%8.3At%10.0At%10.0At%10.0At%10.0f 
\n’  ,cosI,cosII,omega_2); 

%  end  %  end  of  ’  for  atom  =  [12]’ 
fclose(^); 

Czas  =  cputime  -  czas; 

%Cz  =  fix(clock); 


#  3te  3|e  9(c  :|e  9]c :|c  a(e  a(c  3(e  3(c  3^  3)c  9|c  ♦♦  3)e  a|c  9|c  3te  :|c  3K  *******  :|c  sfe  *  )ie  *  3ft  9|c  sK 

J(ea|esfe)fc9ie)fe*3fesfe3|e9ie9fe:fcsfe3f(:)e:|c:ie3)c***  *********************  ************** 

%  ***  hmx_2o.m  **’“  Jolanta  M.  Snowel  July  19,  1995. 

% 

%  This  program  calculates  second  order  contributions  to  the  second  moment  <  w^2  > 
%  from  Unlike  Spins. 

%  Nitrogen  nuclei:  spin  1=1,  asymmetry  parameter  !=  zero, 

%  no  magnetic  field. 

%  LAST  Version  :  sec_mom.j21,  June  21,  1995. 

%  Two  sites  are  considered  : 

%  I  -  N(2)  in  the  positions:  (x,y,z)  &  (-x,-y,-z) 

%  n  -  N(2)  in  the  positions:  (-x,y-l/2,-z-l/2)  &  (x,-y+l/2,z+l/2) 

%  (-x,y-l/2,z-l/2)  =  -(x,-y+l/2,z+l/2) 

%  <'w*'2>  eqs  :  2.16  &  3.27  S.Vega 

clear  all 

czas  =  cputime; 


48 


hmx_b  %  hmx_b.ni  file.  It  loads  data:  unit  cell  parameters,  x,  y,  z,  and  B_y  (beta), 
t  =  3;  %  Translation  limits  for  n,  m,  and  p. 

atom  =  1;  %  1  for  N(2)  ;  2  for  N(3) 

fp  =  fopen(’c:\jolanta\hmx\hmx_2o.dat’,’w’); 

%for  atom  =  [1  2] 
count2  =  0; 

^  ^ittitiiiciic*****************  efg  axes  ************************* 

%  Make  a  table  of  pairs  of  atom  numbers  in  each  row, 

%  numbers  refer  to  the  sequence  in  ’xyz’  matrix  in  hmx_b.m  file 
%  (also  see  efg.m,  efg_uc.m,  &  efg_ucat.m): 

at  =  [ 

6  5  12  %  EFGforN(2) 

8  7  4  3  ];  %  EFGforN(3) 

at  =  at(atom,:); 

sqrt_g  =  sqrt(det(g)); 

k  =  size(at,l);  efg  =  zeros(k*12,3);  RR  =  zeros(k*12,3); 

vectorl  =  zeros(4*k,3);  vector2  =  zeros(4*k,3);  vector3  =  zeros(4*k,3); 

for  i  =  1  :  k 

vecl  =  xyz(at(i,2),:)  -  xyz(at(i,l),:);  %  [EFG]1  =  [O  -  O]  <--  auxUiary 

vec2  =  xyz(at(i,4),:)  -  xyz(at(i,3),:);  %  [EFG]2  =  [N  -  N] 

vecll=  sqrt(vecl*g*vecr); 
vec21=  sqrt(vec2*g*vec2’); 

vecl  =  vecl/vecll; 
vec2  =  vec2/vec21; 

vector2(i,:)  =  vec2; 

^  [EFG]3  =  [O  -  O]  X  [N  -  N]  ♦*♦♦******♦♦ 

vec3(l)  =  sqrt_g*(  vecl(2)*vec2(3)  -  vecl(3)*vec2(2) )  ; 
vec3(2)  =  sqrt_g*(  vecl(3)*vec2(l)  -  vecl(l)*vec2(3) )  ; 
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vec3(3)  =  sqrt_^’''(  vecl(l)*vec2(2)  -  vecl(2)*vec2(l) )  ; 


vector31  =  sqrt(vec3*G*vec3’); 
vector3(i,:)  =  (vec3/vector31)*G  ; 
vec3  =  vec3*G  ; 

^  [EPG]1  =  [EFG]2  X  [EFG]3  ************* 


vectorl(i,l)  =  sqrt_g*(  vec2(2)*vec3(3)  -  vec2(3)*vec3(2) )  ; 
vectorl(i,2)  =  sqrt_g*(  vec2(3)*vec3(l)  -  vec2(l)*vec3(3) )  ; 
vectorl(i,3)  =  sqrt_g*(  vec2(l)*vec3(2)  -  vec2(2)*vec3(l) )  ; 

vectorl(i,:)  =  vectorl(i,:)*G  ; 

vectorll  =  sqrt(vectorl(i,:)*g*vectorl(i,:)’); 

vectorl(i,:)  =  vector  l(i,:)/vectorll; 

end  %  *********’►**  END  ************** 

%  ******  EFG  axes  of  all  N(2)  and  N(3)  in  equivalent  positions  ******* 

vectorl(k+l:2*k,:)  =  [-vectorl(l:k,l)  vectorl(l:k,2)  -vectorl(l;k,3)];  %  2) 

-x,y+l/2,-z+l/2 

vector2(k+l:2*k,:)  =  [-vector2(l:k,l)  vector2(l:k,2)  -vector2(l:k,3)];%  +  0.5 

9 

vector3(k+l:2*k,:)  =  [-vector3(l:k,l)  vector3(l:k,2)  -vector3(l:k,3)];%  +  0.5 


vectorl(2*k+l:3*k,;)  =  -vectorl(l:k,:)  ;  %  (3)  -x,-y,-z 

vector2(2*k+l:3*k,:)  =  -vector2(l:k,;)  ; 

vector3(2*k+l:3*k,;)  =  -vector3(l:k,:)  ; 

vectorl(3*k+l:4*k,:)  = -vectorl(k+l:2*k,:);  %  (4)  x,-y+l/2,z+l/2 

vector2(3*k+l:4*k,:)  =  -vector2(k+l:2*k,:);  %  +  0.5  ; 

vector3(3*k+l:4*k,:)  =  -vector3(k+l:2*k,:);  %  +  0.5  ; 

^  itutcUcitcilt DC 

E  =  [vectorl(l,:)  ;  vector2(l,:)  ;  vector3(l,:)];  invE  =  inv(E)  ;  %  EFG_x 

%  E  =  [vectorl(3,:)  ;  vector2(3,:)  ;  vector3(3,:)];  invE  =  inv(E)  ; 

n=l;  1=3; 


for  i  =  1  :  4*k 
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efg(n:l,:)  =  [vectorl(i,:)  ;  vector2(i,:)  ;  vector3(i,:)]; 

F  =  efg(n:l,:); 

RR(n:l,:)  =  E*inv(F); 
n  =  l+l;  1  =  1  +  3; 
end 

%  ♦*♦*♦*♦♦♦♦♦  end  of  calculations  of  the  efg  tensors 

sfe  #  3(e  aK  ♦  #  4e  #  sK  3|e  #  4c  3|c 

%  ♦♦♦♦♦**♦**♦  &  TRANSFORMATION  MATRIX  [RR] 

4c  #  4e  4c  4e  4e  4(  ife  ♦  9k  4e  4c  4e  )|c  4c  ♦  aK  4c  4e  4e  ♦  4e  ♦  4<  4«  4c 


% 

4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4e4c4c4c4c4c4c4e4c4c4e4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c4c 

4c4c4c4c4c4c4c 

%  ♦♦**♦♦♦♦♦♦  POSITIONS  OF  EQUIVALENT  NUCLEI  {  N(2)  } 

4c  4c  4c  4c  4c  4c  4c  4c  4e  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c  4c 

%  Either  atoms  N(2)  or  N(3)  from  ’xyz’  matrix  (hmx  b.m  file)  are  considered: 
m  =  5  ; 

xyz_d  =  zeros(m*4,3); 

xyz_d(l,:)  =  xyz(atom+l,:);  %  fractional  coordinates  of  N(2)  or  N(3) 
xyz_d(2:5,:)  =  xyz(ll:14,:);  %  fractional  coordinates  of  protons 


%  ***  xyz_d  -  matrix  of  the  coordinates  of  nitrogen  N(2)  atoms 
%  ***  in  the  following  positions: 

%  ♦**  (1)  x,y,z  (2) -x,y-l/2,-z-l/2  (3) -x,-y,-z  (4)  x,-y+l/2,z+l/2  . 

%  ♦**  (1)  &  (3)  -  1st  molecule;  (2)  &  (4)  -  2nd  &  3rd  molecule; 

xyz_d(m+l:2*m,l)  =  -xyz_d(l:m,l)  ;  %  2)  -x,y-l/2,-z-l/2 

xyz_d(m+l:2*m,2)  =  xyz_d(l:m,2)  +  0.5  ; 

xyz_d(m+l:2*m,3)  =  -xyz_d(l:m,3)  +  0.5  ; 


:%xyz_d(m+l:2*m,l)  = -xyz_d(l:m,l)  ;  %  2) -x,y+l/2,-z+l/2 

%xyz_d(m+l:2*m,2)  =  xyz_d(l:m,2)  -  0.5  ; 
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%xyz_d(m+l:2*m,3)  =  -xyz_d(l:m,3)  -  0.5 


xyz_d(2*m+l:3*m,l)  = -xyz_d(l:m,l)  ;  %  (3) -x,-y,-z 

xyz_d(2*m+l:3*ni,2)  =  -xyz_d(l:m,2)  ; 

xyz_d(2*m+ 1:3*111,3)  =  -xyz_d(l:m,3)  ; 

xyz_d(3*m+l:4*m,l)  =  xyz_d(l:m,l)  ;  %  (4)  x,-y+l/2,z+l/2 

xyz_d(3*m+l:4*m,2)  =  -xyz_d(l:m,2)  +  0.5  ; 

xyz_d(3*iiH-l:4*m,3)  =  xyz_d(l:m,3)  +  0.5  ; 


k  =  size(efg,l);  countl  =  0;  count2  =  0; 
i=l; 

for  j  =  l:(k/3) 
vic2(j,:)  =  efg(i,:)*invE; 
i=i+3; 
end 

%  *****♦♦♦♦♦  SECOND  MOMENT  CONTRIBUTION  FROM  UNLIKE  SPINS  ♦*♦♦♦♦***** 
SecMoment  =  zeros(6,l);  SecMom  =  0;  %SecondMom  =  0; 

^  bCtWCCfl  2  atOITlS  ♦***************’t^******* 

k  =  size(xyz_d,l);  countl  =  0  ; 

%  vi  spin  I  (nitrogen)  %  vs  (hydrogen)  %  vt  (hydrogen) 


%for  i  =  1:5:20 
for  i  =  [1  6] 

sm_NH  =  0  ; 
vi  =  xyz_d(i,:); 

for  j  =  [2:5  7:10  12:15  17:20] 

fori  =  [2:5  7:10  12:15  17:20] 


52 


ifa  >  j) 


% 

vs  = 

xyz_d0.0; 

% 

vt  = 

xyz_d(l,:); 

for 

ms 

=  -t;t 

vs(l) 

=  xyz_d(j,l) 

+ 

ms  ; 

for 

ns 

=  -t:t 

vs(2) 

s 

X 

II 

+ 

ns  ; 

for 

ps 

=  -t:t 

vs(3) 

>> 

X 

II 

+ 

ps  ; 

for 

mt 

=  -t:t 

vt(l) 

=  xyz_d(l,l) 

+ 

ms  ; 

for 

nt 

=  -t:t 

vt(2) 

=  xyz_da,2) 

+ 

ns  ; 

for 

Pt 

=  -t:t 

vt(3) 

=  xyz_d(l,3) 

+ 

ps  ; 

countl  =  countl  +1  ; 


v_s  =  vs  -  vi ; 

v_t  =  vt  -  vi  ; 

r_s  =  sqrt(v_s*g*v_s’)  ; 

r_t  =  sqrt(v_t*g*v_t’)  ; 

D_s  =  v_s*invE  ; 

X_s(l)  =  D_s(l)/r_s;  Y_s(l)  =  D_s(2)/r_s;  Z_s(l)  =  D_s(3)/r_s; 

XOS_s(l)  =  1  -  3*D  s(l)*D_s(l)/(r_s*r_s); 

YOS_s(l)  =  1  -  3*D_s(2)*D_s(2)/(r_s*r_s); 

ZOS_s(l)  =  1  -  3*D_s(3)*D_s(3)/(r_s*r_s); 

D_t  =  v_t’'‘invE  ; 

X_t(l)  =  D_t(l)/r_t;  Y_t(l)  =  D_t(2)/r_t;  Z_t(l)  =  D_t(3)/r_t; 

XOS_t(l)  =  1  -  3*D  t(l)*D_t(l)/(r_t*rJ); 

YOS_t(l)  =  1  -  3*D_t(2)*D_t(2)/(r_t*r_t); 

ZOS_t(l)  =  1  -  3*D_t(3)*D_t(3)/(r_t*r_t); 
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%  ***  Second  Order  Contributions  from  Unlike  Spins.  **** 
%  ***  eq.  3.27  *** 


sm  =  (2*ZOS_s*ZOS_s*ZOS_t»ZOS_t  +  ... 

9*(ZOS_s*ZOS_s*(X_t*X_t  +  Y_t*Y_t)  +  ZOS_t*ZOS_t»(X_s*X_s  + 
Y_s*Y_s))  +  ... 

81*(X_s*X_s*X_t*X_t  +  Y_s*Y_s*Y_t*Y_t  +  ... 

X_s*X_s’"Y_t*Y_t  +  Y_s^Y_s*X_t*X_t)  )/r_s^6/r_r6  ; 

sm_NH  =  sm_NH  +  sm  ; 

%  Distance(countl,:)  =  [i  j  1  r_s  r_t  sm]  ; 

%fprintf(fp,’%6.0f  -%3.0f -%3.0At  %8.4f  %10.4f  %22.10f  \n’,[i  j  1  r_s  r_t  sm]); 


end  %  end  ’  for  pt  =  -t:t  ’ 

end  %  end  ’  for  nt  =  -t:t  ’ 

end  %  end  ’  for  mt  =  -t:t  ’ 

end  %  end  ’  for  ps  =  -t:t  ’ 

end  %  end  ’  for  ns  =  -t:t  ’ 

end  %  end  ’  for  ms  =  -t:t  ’ 

end  %  end  ’  if  G  >  j)  ’ 

end  %  end  ’  for  1  =  [2:5  7:10  12:15  17:20]  ’ 

end  %  end  ’  for  j  =  [2:5  7:10  12:15  17:20]  ’ 

count2  =  count2  +1; 

SM_NH(count2)  =  sm_NH; 

end  %  end  ’  for  i  =  1:5:20  ’ 

%  const  =  (1.9338*26.7520*1.05459)"4/(16*pi*pi*1564*1564)*1000000  ; 
%  N(3) 

const  =  (1.9338*26.7520*1.05459)^4/(16*pi*pi*1441*1441)*1000000  ; 
%  N(2) 

SM_NH_av  =  0.5*sum(SM_NH*const);  %  const  = 
(l/4)*gamma_N'^4’''gamma_H'^4’'‘h_bar''6/(E_x  -Ejy)''! 


%fprintf(’%6.0f  -  %3.0At  %3.0f  %3.0f  %3.0f  %12.4f%12.4f%12.4f%20.8f 

\n’, Distance’); 
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^  Wilting  tO  E  filC  ****************** 

Q)rintf(fp,’Data  for  HMX  ;  N(%1.0f)\t  t  =  %1.0f;\t  countl  =  %12.0f  \t  Dell 
\n’,atom+ l,t,countl); 

fprmtf(fp,’\n\n  SM_NH  =\t%20.6f  %20.6f  %20.6f  %20.6f  \n’,  SM_NH); 
fprintf(fp,’\n\n  SM_NH_av  =%20.6f  sqrt(SM_NH_av)  =  %20.6f  \n\n  \n’,  [SM_NH_av 
sqrt(SM_NH_av)]); 

%end 

fclose(fp); 

SM_NH’ 

SM_NH_av 

Czas  =  q)utime  -  czas 
%Cz  =  fix(clock); 

%save  last_s.mat 

♦*♦*♦♦♦♦♦**♦♦***♦♦♦»♦♦♦♦*♦♦*♦*♦**»»♦**♦**♦********♦**************** 

♦♦♦♦♦#»#*#»♦♦♦♦»♦»*♦♦♦♦***♦»♦♦♦*»**♦»♦♦»*»**♦*♦***********♦♦***♦♦** 
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